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Abstract 

This paper introduces a new way of generalizing Hilbert's two-dimensional space-filling curve to 
arbitrary dimensions. The new curves, called harmornous Hilbert curves, have the unique property 
that for any d! < d, the d-dimensional curve is compatible with the d'-dimensional curve with 
respect to the order in which the curves visit the points of any d'-dimensional axis-parallel space 
that contains the origin. Similar generalizations to arbitrary dimensions are described for several 
variants of Peano's curve (the original Peano curve, the coil curve, the half-coil curve, and the 
Meurthe curve). The d-dimensional harmonious Hilbert curves and the Meurthe curves have neutral 
orientation: as compared to the curve as a whole, arbitrary pieces of the curve have each of d! possible 
rotations with equal probability. Thus one could say these curves are 'statistically invariant' under 
rotation — unlike the Peano curves, the coil curves, the half-coil curves, and the familiar generalization 
of Hilbert curves by Butz and Moore. 

In addition, prompted by an application in the construction of R-trees, this paper shows how to 
construct a 2d-dimensional generalized Hilbert or Peano curve that traverses the points of a certain 
d-dimensional diagonally placed subspace in the order of a given d-dimensional generalized Hilbert 
or Peano curve. 

Pseudocode is provided for comparison operators based on the curves presented in this paper. 

1 Introduction 

Space-filling curves A space-filling curve in d dimensions is a continuous, surjective mapping from 
M to M''. In the late 19th century Peano [15] described such mappings for d = 2 and d = 3. Since 
then, quite a number of space-filling curves have appeared in the literature, and space-filling curves have 
been applied in diverse areas such as spatial databases, load balancing in parallel computing, improving 
cache utilization in computations on large matrices, finite element methods, image compression, and 
combinatorial optimization 2 . 

Space-filling curves are usually recursive constructions that map the unit interval [0, 1] to a subset 
of M.'^ that has measure larger than zero. In the case of Peano's two-dimensional curve, the subset of 

that is 'filled' is the unit square. Peano's curve is based on subdividing this unit square into a 
grid of 3 X 3 square cells, and simultaneously subdividing the unit interval into nine subintervals. Each 
subinterval is then matched to a cell; thus Peano's curve traverses the cells one by one in a particular 
order. The mapping from unit interval to unit square is refined by applying the procedure recursively to 
each subinterval-cell pair, so that within each cell, the curve makes a similar traversal (see Figure [ija)). 
The result is a fully-specified mapping from the unit interval to the unit square. A mapping from M to 

could then be constructed by inverting the recursion, recursively considering the unit interval and 
the unit square as a subinterval and a cell of a larger interval and a larger square. A higher-dimensional 
version would be based on subdividing a d-dimensional hypercube into a grid of S"* hypercubic cells. 

*This manuscript has been put on arXiv.org for early dissemination of the results. I have not found the time (yet) 
to prepare this manuscript (or parts of it) for submission to be published in a journal or a conference. I welcome any 
corrections and other comments which could help me to improve this manuscript and/or to disseminate the results further: 
please mail me at cs.herman@haverkort.net. 

iDept. of Computer Science, Eindhoven University of Technology, the Netherlands, cs.herman@haverkort.net 
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Figure 1: (a) A sketch of Peano's space-filling curve, (b) A sketch of Hilbert's space-filling curve. 



Peano's curve has been the curve of choice for certain applications. However, in many other appli- 
cations preference is given to curves based on subdividing squares (or hypercubes) into only 2'* squares 
(or hypercubes) . The cell in which a given point p lies can then be determined by inspecting the binary 
representations of the coordinates of p bit by bit, and no divisions by three need to be computed. In 
response to Peano's publication, Hilbert [TU] described a two-dimensional space-filling curve based on 
subdividing a square into four squares (Figure [ijb)). However, he did not describe how to do something 
similar in higher dimensions. It is now well-established that it can be done; in fact there are many ways 
to define a three-dimensional curve based on subdividing a cube into eight octants [HIS]. 



The goal: extradimensional space-filling curves In this paper we consider different ways of gen- 
eralizing Hilbert curves, Peano curves and similar curves to higher dimensions. Let P be a set of points 
on a d'-dimensional facet U' of a d-dimensional unit hypercube U, such that the origin is a vertex of U 
and U' . In this paper, we study generalizations of space- filling curves that have the following property: 
regardless of our choice of P and U', the order in which the points of P are visited by the d-dimensional 
curve that fills U, is the same as the order in which these points are visited by the d'-dimensional curve 
that fills U'. 

To get some intuition for what this means, consider Figure [2j Figure [2ja) shows the order in which 
the two-dimensional Peano curve traverses the nine subsquares of the unit square, and Figure [2|^b) shows 
the order in which a three-dimensional Peano curve traverses the 27 subcubes of the unit cube. If we 
consider only the subcubes that appear on the front, left, or bottom face of the unit cube, then we see 
that on each face, these subcubes are visited in the order of the two-dimensional Peano curve. FigurejSja) 
shows the order in which the two-dimensional Hilbert curve traverses the four subsquares of the unit 
square, and Figure |3][b) shows the order in which the three-dimensional generalization by Moore [H] 
(citing Butz [S]) traverses the eight subcubes of the unit cube. Again, we see that on each face, the 
subcubes are visited in the order of the two-dimensional curve. However, from these observations one 
cannot conclude that all points of the front, left or bottom face are visited in the order of the two- 
dimensional curve. Figure |4] shows the next level of refinement of the two-dimensional Hilbert curve and 
the three-dimensional Butz-Moore curve. The order in which the subcubes of the front, left or bottom 
face are visited no longer corresponds to the two-dimensional curve. Our goal in this paper is to identify 
generalizations of the Hilbert curve and several variants of the Peano curve, in which the correspondence 
is maintained for any level of refinement of the order in which squares and cubes are visited. 

To make this statement more precise, we should first define the order in which points appear along the 
curve more precisely. Let / be a space-filling curve, or more precisely, a continuous surjective function 
from [0, 1] to a unit hypercube U. Let a valid ordering / of / be any bijection from the unit hypercube 
U to [0, 1], such that for any point p £ P we have f{f{p)) = P- 

Let /i : {0,...,d' — 1} {0,...,d— 1}, d > d' he a, strictly increasing injection, that is, for any 
i,j G {0, — 1} with i < j we have fi(i) < ^{j). This function will be used as a function that 
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Figure 2: The Peano curve in two and three dimensions, (a) Ordering of nine squares in 2D. 
(b) Ordering of 27 cubes in 3D, with the induced orderings on the front, left, and bottom side. 
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Figure 3: The Hilbert curve and the three-dimensional version of Butz-Moore. (a) Ordering of four 
squares in 2D. (b) Ordering of eight cubes in 3D, with the induced orderings on the front, left, and 
bottom side. 



selects d' dimensions out of d, in the following way: given a rf'-dimensional point p = (po, ...,Pd'-i), we 
define lift^(p) as the d-dimensional point ((Jq, ...,<7(i-i) such that qj = pi when j = fi{i), and qj = when 
there is no i such that j = In other words, lift^(p) is obtained from p by inserting zero coordinates 
at the positions not selected by /x. For a set of points P, let lift^(P) be the set (J^^p lift^(p). 

Definition 1 A d-dimensional space-filling curve f is extradimensional to a d' -dimensional space-filling 
curve f with co-domain U', if the following holds for any valid ordering f of f and any choice of /x; 
there is a valid ordering f of f such that for any pair of points a,b G U' we have f'{a) < f'{b) if and 

onlyifJ{m^,ia))<J{mt^{b)). 

Definition 2 An (infinite) set F of space-filling curves is interdimensionally consistent, if F contains 
a unique d-dimensional space-filling curve fd for any integer d> 1, and fj G F is extradimensional to 
fiGF whenever j > i. 
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Figure 4: The Hilbert curve and the three-dimensional version of Butz-Moore. (a) Ordering of 16 squares 
in 2D. (b) Ordering of 64 cubes in 3D, with the induced orderings on the front, left, and bottom side. 

Inverting the terminology, we could say that /^j/ is intradimensional to fd if and only if is ex- 
tradimensional to fd', and each curve of an interdimensionally consistent set could be regarded as an 
intradimensional curve of an infinitely-dimensionaj^ 'curve' foo. 

In the definition of extradimensionality, a d-dimensional point is constructed from a c?'-dimensional 
point by inserting zero coordinates at the positions not selected by /x. We will also consider first-half- 
extradimensionality, where a 2(i-dimensional point is constructed by appending d zeros to the coordinates 
of a d-dimensional point, and diagonal- extradimensionality, where a 2(i-dimensional point is constructed 
by concatenating two copies of the coordinate sequence of a d-dimensional point. In our previous work 
we found that first-half- and diagonal-extradimensional space-filling curves can be useful in improving 
the performance and robustness of certain data structures based on space-filling curves |9j . 

The results This work provides descriptions of interdimensionally consistent generalizations of Hilbert's 
curve and several variants of Peano's curve (coil, half-coil, and Meurthe) that caught my attention in 
previous work [71 |S] because of their excellent locality-preserving properties. Furthermore, this paper 
provides a general technique to derive first-half- and diagonal-extradimensional curves from a given gen- 
eralized Hilbert or Peano curve. This enables the extension of our previous work on R-trees for rectangles 
to higher dimensions [9] . 

It turns out that the new generalized Hilbert curves and the generalized Meurthe curves share another 
interesting property: they are, in some sense, statistically invariant under rotation. For an intuitive 
understanding of what this means, one could look at Figure [l] At the highest level, the Peano curve 
shows more vertical steps than horizontal steps, and in recursion, it stays that way. Thus the curve has 
a clear orientation, which can be recognized immediately even if one sees only a small piece of the curve. 
This could be a disadvantage if the recognizable orientation somehow interacts negatively with the data 
in an application of the curve. On the other hand, at the smallest depicted level of detail of the Hilbert 
curve, all possible orientations of the curve occur approximately equally often. Therefore, from the looks 
of a small piece of the Hilbert curve one could not tell if the curve as a whole has been rotated or not: we 
say the curve has neutral orientation. The harmonious Hilbert curves and the Meurthe curves maintain 
this property also in higher dimensions — unlike the Peano curves, coil curves, half-coil curves and the 
Butz-Moore generalization of the Hilbert curve. 

have not explored where the idea of infinitely-dimensional curves would lead us. I believe they might constitute a 
new class of "monster curves" (to use Gardner's words |4J). 
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How to read this paper Below, in Section [2j we fill first see how exactly we can define and describe 
space-filling curves like Hilbert's and Peano's. We will distinguish 2-regular curves, which are based 
on subdividing d-dimensional hypercubes into 2"* subcubes, and 3-regular curves, which are based on 
subdividing d-dimensional hypercubes into 3"^ subcubes. In Section [s] we see some basic results that will 
come in useful in proving the claimed properties of the curves presented in the sections that follow. In 
Section [4] we will present interdimensionally consistent generalizations of several 3-regular curves, and 
prove that among those, the generalized Meurthe curves have neutral orientation. We will also give 
pseudocode for comparison operators that decide which of any two given points comes first on the curve. 
The pseudocode may not be a recipe for highly optimized code in any particular programming language, 
but it may give some insight in how one may produce such code and, in addition, it reveals properties 
of the curves that would be difficult to see otherwise. In Section [5] we learn about the Butz-Moore 
generalization of Hilbert curves and about our harmonious Hilbert curves in detail. Again, pseudocode 
is included. Section [6] defines the concepts of first-half- and diagonal-extradimensionality in detail, and 
explains how, given any 6-regular curve / (for b £ {2,3}), one can construct a ^-regular curve g that is 
first-half- and/or diagonal-extradimensional to /. The paper concludes with Section [T] with a comparison 
between the results from the different sections of this paper and with open questions for possible further 
research. 

In principle. Sections [4j [5] and [6] are independent from each other. The reader who wants to learn 
quickly how to get interdimensionally consistent, first-half-extradimensional or diagonal-extradimensional 
curves, may read Section [2] and then skip to any of the Sections |4] [5] or [6] right away. However, the reader 
who wants to understand the proofs in these sections, is advised to read Sections [2] [3] |4] and [5] in order, 
because the proofs in these sections build up in complexity. 



2 Describing space-filling curves 

2.1 Graphical descriptions and their interpretation 

In general, we can define an order (scanning order) of points in d-dimensional space as follows. We give a 
set of rules, each of which specifies (i) how to subdivide a region in d-dimensional space into subregions; 
(ii) what is the order of these subregions; and (iii) for each subregion, which rule is to be applied to 
establish the order within that subregion. We also specify a unit region, and we indicate what rule is 
used to subdivide and order it. 

Figure [5 'a) shows an example (the scanning order depicted corresponds to a section of the f3Q-cuTve 
by Wierum TB]). Each rule is identified by a letter, and pictured by showing a unit square, its subdivision 
into quadrants, the scanning order of the quadrants (by a directed curve along the outer vertices of the 
quadrants), and the rules applied to the quadrants (by letters). Within each quadrant, there is a scaled 
copy of the base pattern: the curve along the outer vertices of the quadrants. Variations of the rules 
that consist of simply rotating or mirroring the order within quadrants are indicated by rotating or 
mirroring the pattern. Variations of rules that consist of reversing the order of cells within a quadrant 
are indicated by reversing the direction of the pattern and by an overscore above the letter identifying 
the rulfj^ Figure [sjjb) shows how applying the rules recursively, starting with rule R, results in an order 
of the subquadrants within each quadrant. Combined with the order of the quadrants, this gives an order 
of all squares in a 4 x 4 grid, which is sketched by the curve shown in Figure [5](c). From the orientations 
of the base patterns in Figure [sjb) and (c), one already gets a hint of the order of all squares in a 8 x 8 
grid. By expanding the recursion further, one may order the cells of an arbitrarily fine grid. 

The polygonal curve that connects the centres of these cells in order thus forms an arbitrarily fine 
approximation of the space-filling curve, which is more precisely defined as follows. At any level of 
recursion, let the grid consist of the cells Ci,...,C„, indexed from 1 to n according to the scanning 



^Reversals are among the standard transformations applied in the definitions of space-filling curves in the literature — see, 
for example, the Gosper flowsnake [4|, arguably one of the prototypical space-filling curves along with the curves by Peano, 
Hilbert, Lebesgue [12J and Sierpinski [2]. In fact, reversals are so standard, that they are often not described explicitly, 
and it is left to the reader to infer the reversals from how the pieces of the curve are put together. I prefer to be explicit 
about reversals. To indicate a reversal, it is necessary to reverse the direction of the pattern and to include the overscore. 
Figure[5ja) illustrates this. Rule J is asymmetric because different rules are applied to the first and the last quadrant. It is 
therefore important to distinguish between refiected, unreversed applications of J, as in the lower left quadrant of rule R, 
and unreflected, reversed applications of J, as in the lower right quadrant. Since the base pattern of J is symmetric, it is 
only the overscore that indicates which of these alternatives is meant. 
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Figure 5: (a) Definition of a scanning order corresponding to a variant of the pQ-cnrve. (b) Applying 
the rules recursively within R. (c) The resulting order of the 4 x 4-grid. (d) Order of the 16 x 16-grid. 
On each level, the last cell of each quadrant shares a vertex (and even an edge) with the first cell of the 
next quadrant. 

order, let vol[l,i] = J2]=i^'^K^i) be the total volume of cells Ci up to and including Q, and define 
fraction [ 1, i] = vol[l, i]/vol[l, n]. Then the curve / maps the interval [fraction[l, i — 1], fraction[l, j]] to 
Ci, for alH G {1, n}. The image f{q) of any point or interval g in [0, 1] can now be obtained by taking 
the limit of the image of the intervals containing q as the level of recursion goes to infinity. Thus the 
definition of a scanning order is also a definition of a space-filling curve. 

In the definitions of space-filling curves, the letters indicating the rules may be omitted if there 
are no reversals and no two rules have the same base pattern modulo rotation and refiection — see, for 
example. Figure [6j Figures [7} [8} [9| and [T0| and Inset [T] show further examples of definitions of two- and 
three-dimensional space-filling curves. 

Entrance and exit gates When the scanning order is refined to an infinitely fine grid, the first and 
the last cell visited shrink to points. These points are the images of and 1 of the space-filling curve 
defined by the scanning order, and they typically lie on the boundary of the unit regiorj^ We will call 
/(O) the entrance gate and we will call /(I) the exit gate of the curve; in Figure [s] their locations are 
indicated by the white dots and the black dots, respectively. Note that these dots are only shown for 
clarity: the locations of the gates are a consequence of the definition of the scanning order, rather than 
a part of the definition. 

Using space-filling curves For many practical applications of space-filling curves it is not necessary 
to actually draw a curve. It is often enough to be able to decide for any two given points p and q, which 
of these appears first on the curve — that is, which point comes first in the scanning order. This can be 
decided by expanding the recursion only to the smallest depth at which p and q lie in different quadrants. 
A technical problem with this is that, down from some depth of recursion, a given point p may always 

^It is not automatic that these points lie on the boundary, but space-fiUing curves are usually designed such that /(O) 
and /(I) do lie on the boundary — otherwise the curve cannot be continuous from one subregion into another subregion if 
each subregion is filled with a scaled-down copy of the curve. 
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lie on the boundary of two or more quadrants. This may create ambiguity about which of two given 
points comes first. This ambiguity can be resolved by using a consistent tie-breaking rule. For example, 
in the pseudocode in this paper, we will always assign each point p to the quadrant that lies to its upper 
right — that is, for a space-filling curve / we define f{p) as the common limit of the values t such that 
f{t) = p -t- (e, e) as e > approaches zero. 

2.2 Classes of space-filling curves considered in this paper 

There are many different possible generalizations of Peano's and Hilbert's curves to three and more 
dimensions. To be able describe the classes of curves we will be studying in this paper precisely, I will 
now introduce some terminology. 

All space-filling curves studied in this paper are defined by a rule system as described above. 

We will call a curve b-regular if the unit region of each rule is a d-dimensional unit hypercube, and 
each rule subdivides this hypercube into b'^ smaller hypercubes of equal size. We adopt the convention 
that the unit hypercube is an axis-parallel hypercube that spans the interval [0, 1] in each dimension; 
thus, each of the smaller hypercubes has width 1/6. 

We will call a curve a mono-curve if the defining rule system contains only a single rule. 

We will call a curve order-preserving if it is a mono-curve and the transformations of the single 
definining rule within the subregions are restricted to rotation and refiection — excluding reversal of the 
order in which the subregions are traversed. 

We will call a scanning order and the space-filling curve defined by it vertex-continuous if it has the 
following property: whenever the scanning order visits a set of hypercubes Ci, C„, in that order, and 
we refine the scanning order in each of these hypercubes recursively to depth £ so that we obtain an 
order on n ■ b^'^ cells, for any i > 1, then the last cell within Cj_i shares at least one point with the first 
cell within Ci, for all 1 < i < n (see, for example. Figures [l] and [5jd)). Note that if a scanning order is 
not vertex-continuous, then it does not directly define a continuous mapping from the unit interval to 
a higher-dimensional space, and thus, the scanning order does not actually define a space-filling citrwe|^ 
One of our main tasks in the following sections of this paper, will always be to prove that whatever is 
presented as a space-filling curve, is indeed vertex-continuous. 

Inset [T] shows some of the variety of vertex-continuous 2-regular mono-curves that exist. 

2.3 Numerical descriptions 

Obviously, the graphical way of defining space-filling curves described above does not extend well to 
curves in more than three dimensions. Therefore, I will now introduce a concise numerical way of 
defining space-filling curves, tailored to 6-regular order-preserving curves. 

In the numerical notation, we will use d-dimensional coordinate vectors and d-digit numbers. We use 
subscripts from to d — 1 to index the individual coordinates or digits, respectively. If a; is a d-digit 
number, xq is the most significant and Xd~i is the least significant digit of x. By (a)d we denote the 
number that consists of d times the digit a. 

The rank of each subregion is a base-fe number from {0, (b — l)d} that indicates the position of 
the subregion in the sorted order of subregions according to the rule that defines the curve: the first 
subregion has rank and the last subregion has rank (b — l)d- By S(r) we denote the subregion with 
rank r. For ease of notation we write f for r — 1, and we write f for r -\- 1. 

The location of S{r) is denoted by c(r), which is a d-digit base-6 number indicating that the subregion's 
corner that is closest to the origin is located at i{co{r), Cd-i{r)) and the subregion's corner that is 
farthest from the origin is located at i{co{r) -t- 1, Cd-i{r) + 1). 

The curve within each subregion S{r) can be obtained from the curve in the unit hypercube as a 
whole by first rotating it, then reflecting it, then scaling it with a factor 1/6 in all dimensions, and then 
translating it into place. I will now describe how rotations, reflections and translations are specified. The 
rotation is specified by a permutation a(r), which is a d-digit base-d number of which the digits form a 
permutation of {0, ...,d — 1}. This permutation specifies which axis ai{r) of the curve through the unit 

work-around would be to insert infinitely many connecting line segments to bridge the gaps. In effect, this is 
how Lebesgue makes a space-filling curve out of the non-vertex-continuous scanning order that is nowadays known as 
Z-order [12| . However, such tricks come at the expense of decreased performance in applications 8 , and therefore they are 
not considered in this paper. 
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Inset 1 What does it take to be a Hilbert curve? 

Hilbert curves are often associated with 2-regular, order-preserving, symmetric mono-curves, that visit the subregions of the 
hypercube in the order of the binary reflected Gray code (see Section |2.5| , and have the entrance and exit gate located at 
the two endpoints of an edge of the hypercube. Indeed, Hilbert's original two-dimensional curve has all of these properties, 
but one may argue that they do not define the curve. There is only one vertex-continuous 2-regular mono-curve, and this 
happens to be the order-preserving, symmetric, reflected-Gray-code-based, vertex-gated Hilbert curve. 

In three-dimensional space, there is not a single vertex-continuous 2-regular mono-curve, but there are 10 694 807 of them — 
not counting rotated, reflected and/or reversed copies of the same curve [6]. We will call all of these curves mono- Hilbert 
curves. Only 48 of them have all of the non-defining properties mentioned above. The 48 order-preserving, symmetric, 
reflected-Gray-code-based, vertex-gated curves include those depicted in Figures (a), (b) and (c) below. Figure (a) shows 
the curve implemented by Moore 1141 . based on the work of Butz. Curve (b) is the only three-dimensional mono-Hilbert 
curve that is extradimensional to the two-dimensional Hilbert curve; it is the three-dimensional harmonious Hilbert curve 
as deflned in Section [s] Curve (c) has a very regular structure: it is not only composed of eight similar octants and of 
two symmetric halves, but also of four quarters (of two octants each) which arc reflected, rotated and/or reversed copies 
of each other. 

The remaining 10 694 759 mono-Hilbert curves include curves that follow an alternative Gray code (e.g. Figure (d)) and 
curves that are not order-preserving (e.g. Figure (e)). There is even one curve that has its entrance and exit gates in 
the interiors of two faces of the cube (Figure (f)): in fact, this is the three-dimensional mono-Hilbert curve that has the 
best locality-preserving properties if measured according to the Euclidean dilation measure [6]. The curves in Figures (e) 
and (f) are the only two mono-Hilbert curves such that on any level recursion, the distance between the centroids of the 
predecessor and the successor of a subregion is only \/2 times the width of the subregion — the centroids of three consecutive 
subregions are never coUinear. There are also curves with 'diagonal connections', where successive subregions only share 
an edge or a vertex: for example, see Curve (g) (with a similar regular structure as Curve (c)) and Curve (h) (rotationally 
symmetric in the left-right axis, it may be more regular than it seems). 




Some three-dimensional vertex-continuous 2-regular mono-curves, with the names I gave them [6]. (a) A26.00.00, the 
three-dimensional Hilbert curve by Butz and Moore; (b) A26.00.db, the three-dimensional harmonious Hilbert curve; 
(c) A16. 00.3c; (d) A6.00.ff; (e) A26.2b.b3; (f) F; (g) B15.00.c3; (h) B138.00.99 (not published before). 



hypercube is rotated onto axis i within subregion r. By a{r) we denote the inverse of this permutation, 
that is, ai{r) = j if and only if aj{r) = i. The reflections are specified by a d-digit binary number m(r), 
where mi{r) = 1 if and only if after rotation, the curve within subregion r should be reflected in a plane 
orthogonal to axis i. The required translation vector o(r) is implied by the rotation, the reflections and 
the location of the subregion, and does not need to be specified separately. 

More precisely, a(r) and m(r) define a transformation matrix M[r) with rows and columns numbered 
from to d — 1 and the following entries: 

-^(^)jj = 1 if i = 0'i{'>') and mi{r) = 0; 
M{r)ij -1 if j = aj(r) and TOi(r) = 1; 
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Figure 6: Hilbert's curve (left) and Peano's curve (right) defined graphically and by means of a table. 



By M(r) we denote the inverse of this matrix: 

M{r)ji ~ 1 if i = aj{r) and niiir) = 0; 
M(r)ji = —1 if i = aj{r) and mi{r) = 1; 
M(r)j, = iii^a,{r). 

The transformation T(r) : p — ^M(r)p + o(r) now maps each point p of the curve through the unit 
hypercube to a point in subregion S{r). Recall that the space-filling curve is a function / from [0,1] 
to the unit hypercube [0, 1]''. Considering the subregion with rank r, the function / maps the do- 
main [r / b'^ , {r + 1) / b'^] to S{r). The transformation T{r) describes the relation between, on one hand, 
/ restricted to the domain [r/b'^,{r + l)/b'^], and on the other hand, / on the domain [0,1]: for 
X e [r/b^, [r + l)/6'*] we have f{x) — T{r){f{sr{x))), where Sr{x) — b'^x — r is the function that scales 
up the domain \r/b'^,{r-\- l)/b'^] to [0, 1]. Conversely, there is a valid ordering / of / such that for any 
point p in the interior of subregion r we have: f{p) = (/(T(r)(p))) , where T(r)(p) = M{r)b{p — o(r)) 
and s7(a;) = (x + r)/b'^. 

A curve can now be described by a table that lists, for each subregion S'(r), the location c(r), the 
permutation a(r), the reflections m(r), and, redundantly (if desired for clarity), the entrance and/or exit 
gates of each subregion. Figure |6] shows the definitions of the two-dimensional curves by Hilbert and 
Peano in graphical notation and in table format. 



2.4 Implementing a comparison operator 

We will know see how we can easily implement a comparison operator for an order-preserving 5-regular 
curve. It is not necessary to read this section to be able to understand the descriptions and proofs in the 
rest of this paper, and it is not necessary to read this section to be able to use the pseudocode presented 
in the rest of this paper. However, it may help in understanding why the pseudocode in Sections |4j|6] 
works. 

Our goal in this section is the following: given two points p, q in [0, 1]"*, whose coordinates are given as 
base-6 numbers, decide whether /(p) < /(9). As an ordering / corresponding to / we choose the ordering 
that assigns each point p to the subregion that lies above/behind/to the right of it in all dimensions — that 
is, we define f{p) as the common limit of the values t such that f{t) — p + {e, e) as e > approaches 
zero. This defines / only for points in [0, l)'*; it is not defined for points of which one or more coordinates 
are equal to 1, and our pseudocode will not deal with such points. One could define / also for points with 
one or more coordinates equal to 1, and one could easily adapt the algorithms presented in this paper 
to handle such coordinates, by representing the number 1 as 0.111... (if 6 = 2) or 0.222... (if b = 3). 



The idea of our basic algorithm (Algorithm 2.1) is as follows. The algorithm takes two points p 



and g, each given as a set of d coordinates (p[0], ...,p[d — 1]) and {q[Q], ...,q[d— 1]), respectively. Assume 
that each coordinate is in [0, 1) and it is given as the string of digits that constitute the fractional part 
of its representation as a base-6 number. The algorithm first removes the first digit of the fractional 
part of each coordinate. This is done with a function ExTRACTFiRSTDiGiT(a:), which removes the first 
digit from the digit string x and returns the removed digit — if ExtractFirstDigit is called on an 
empty string, it simply returns a zero. Now let Cp and Cg be the concatenated removed digits of p and q, 
respectively. Determine Vp and rq such that c(rp) — Cp and c{rq) = Cq] thus Vp and rq are the ranks of 
the regions containing p and q, respectively. If rp < Vq, return true; if > Tq, return false; otherwise, we 
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apply the transformation x — >■ M(rp)x to the points p and q that remained after removing the first digits 



of their coordinates, and recurse on the resulting points p and q. In fact, in Algorithm 2.1 we do not 
apply the transformation right away, but instead we maintain a transformation matrix T which is the 
composition of all transformations M{rp) that should have been applied at all higher levels of recursion. 
We apply this composed transformation T only to the digits extracted, on Line [5] 

Algorithm 2.1: General comparison operator for a 6-regular curve: the main idea 

Input: Points p — {p[0], — 1]) and q = {q[0], q[d — 1]), with coordinates in [0, 1), given by 

the digits of the fractional part in base-5 representation 
Output: true if f{p) < f{q), false otherwise 

1 T 4— identity matrix of size d 

2 repeat 
for z -s— to c? — 1 do 

|_ Cp[i] EXTRACTFlRSTDlGIT(p[i]); Cq[i] ^ EXTRACTFlRSTDlGIT(q[i]) 



Cp i T • Cp , Cq i T • Cq 

add 6 — 1 to all elements of Cp and Cg that are negative 

compute rp,rq such that c{rp) = Cp and c(rg) = Cq 

if Tp < Tq then return true else if Vp > rq then return false 



// Cp = Cq and Vp = 
T ^ M{rp) ■ T 

10 until all remaining digit strings p[Q]^ ...^pld — 1] and g[0], ^[(i — 1] are empty 

11 return false / / p and q are equal 



The correctness of Algorithm 2.1 is based on the following observations. To compare f{p) 
{f{T{r){p))) and f{q) = (/(T(r)(g))) we do not have to apply s^ if both points He in the same 
region S{r): we can simply compare f{T{r){p)) to /(r(r)(q)). Algorithm 



2.1 



does this as follows. Re- 
moving the first digits of the coordinates implements the transformation x — >■ b{x — o'(r)), where o'(r) is 
the corner of S{r) closest to the origin. Thus, assuming o{r) = o'{r), subsequently applying x — M{r)x 
implements T(r). If there are reflections, that is, negative entries in M(r), then the situation is slightly 
more complicated. The translation vector o(r) should be the coordinates of the corner of region S{r) 
that maps to the origin under the transformation T(r). Mapping o'(r) to the origin instead, results in 
a translation of the unit hypercube T{r){S{r)) over a distance -1 in each of the reflected coordinates. 
Line |6] corrects this by adding 5 — 1 to each extracted digit of a negative coordinate — which amounts to 
adding 0.111... = 1 (if 6 = 2) or 0.222... = 1 (if 6 = 3). 

Algorithm |2.2| shows a version of Algorithm |2.1| which is more explicit about how to maintain and 
apply the transformation matrix T. Note that T is always a matrix with exactly one non-zero entry in 
each row, and this entry is either 1 or —1. Therefore, instead of maintaining a complete matrix, we can 
simply maintain where the non-zero entries are and whether they are 1 or —1. 



2.5 Reflected Gray codes 

Many variants of Peano's and Hilbert's curves are based on Gray codes: sequences of base-6 numbers of 
d digits such that each number from {0, b'^ — 1} occurs exactly once and consecutive numbers differ in 
only one digit. More precisely, we require that the difference in this digit is always exactly one. When the 
numbers of such a Gray code are interpreted as the locations c{r) of the regions S{r) for r = 0, 6"^ — 1, 
the result is a space-filling curve in which each visited region S{r) shares a (d— l)-dimensional face with 
its predecessor S{f) (recall that we use f to denote r — 1); for examples, see Figures (a)-(f) in Inset [ij 
Thus, each region fits snugly together with the previous region. This tends to result in space-filling 
curves with good locality-preserving properties, although by itself, it is not a guarantee [6]. 

In many of the definitions of space-filling curves that follow, we will use a specific type of Gray codes 
called reflected Gray codes. By 5rgc'* we denote the reflected Gray code with base b. This code is a se- 
quence of b'^ base-6 numbers, containing each of the 6'' different base-6 numbers exactly once. The code is 
defined recursively as follows. Let a®6RGC'^ denote the sequence that is obtained by prefixing each num- 
ber of bRGC'^ with the digit a, let rev(6RGC'') denote 6rgc'^ in reverse order, and let revifodd(a, 6rgc'*) be 
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Algorithm 2.2: General comparison operator for a 6-regular curve: more detailed algorithm 

Input: Points p — (j3[0], ...,p[d — 1]) and q = ((/[O], q\d — 1]), with coordinates in [0, 1), given by 

the digits of the fractional part in base-5 representation 
Output: true if /(p) < /(g), false otherwise 

1 for i -i— to d — 1 do reflected[i] -s— false // maintains which columns of T have minus 

2 for i -i— to d — 1 do aa:is[i] -s— i // maintains which entry in each row of T is non-zero 

3 initialize an array altAxis[0..d — 1] // for a new value of axis under construction 

4 repeat 
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for z ^ to d — 1 do 

Cp[i] <— EXTRACTFlRSTDlGlT(p[aa;is[i]]); 
if reflected[axis[i]] then 



Ca I 



^ ExTRACTFlRSTDlGlT((j[aa;is[i]]) 



[_ Cp[i\ ^ (6 - 1 - Cp[i\); Cg[i\ ^ {b - 



1 



Co I 



compute TpjTq such that c{rp) = Cp and c{rq) = Cq 



if Tp < Tq then return true else if r 



p > Tq then return false 



// Cp = Cq and Vp = r, 

11 for z -s— to d — 1 do 

12 j^a,{rp) 

13 altAxis[i] ^ axis[i] 

14 if mj{rp) = 1 then reflected[altAxis[i\ 

swap axis and alt Axis 

16 until all remaining digit strings p[0], ...,p[d — 

17 return false 



^ -^reflected[altAxis[i 
1] and q[0], (/[d 



1] are empty 

1 1 p and are equal 



6rgc'* if a is even, and rev(6RGC'^) if a is odd. Then 6rgc^ is the sequence (0, 1, 6— 1), and 6rgc'*+-^, 
for d > 1, consists of the concatenation of the sequences a ® revifodd(a, 6rgc'^) for a ~ 0, 1, 6 — 1 in 
order. For example, aRGC^ = (0, 1, 2) and Srgc^ = (00, 01, 02, 12, 11, 10, 20, 21, 22). We can index foRGC'^ 
by d-digit base-5 numbers. More precisely: we can consider ^RGC*^ to be a bijection between the d-digit 
base-6 numbers and the d-digit base-5 numbers, where 6RGC''(r) is the element r (counting from zero) of 
the sequence 6RGC''. Observe that 6rgc''(0) = 0. 

The reader may verify that by construction, two consecutive elements 6RGC''(f) and 6RGC'^(r) differ 
in only one digit, and the difference in value of this digit is one. 

Furthermore, &RGC'^(r) can be obtained from r by going through the digits of r in order from tq to 
r^-x. whenever a digit r,j is encountered such that is odd, we 'reflect' all following digits, that is, we 
replace by 6 — 1 — rj for all j > i. To understand why this works, consider that after having seen r^, we 
have fixed the first i + l digits of 6RGC''(r); if bRGcf{r) is even, then the remaining digits must constitute 
the element of 6rgc'^^'~^ indexed by the base-6 number rj+i...r(i_i, whereas if bRGcf{r) is odd, then the 
remaining digits must constitute the element of rev(6RGC'*^'^^) indexed by the base-6 number ri+i...rd-i, 
or equivalently, the element of bRGG'^~^~^ indexed by r^_|_i...r^_^, where = 6 — 1 — rj 



Conversely, we compute r from 6RGC'^(r) wit h Al gorithm 2.3 Fo r a 6-regular curve of which c(r) is 



defined as 6rgc (r), we can integrate Algorithm 2.3 into Algorithm 2.2 (which computes which of two 



points comes first along a space-filling curve) and obtain Algorithm 2.4 



Algorithm 2.3: Gray-decoding algorithm for 6RGC'' 



Input: An element of 6rgc'^, given as an array c = (c[0], c[d — 1]) of base-6 digits 
Output: rank s.t. bRGC'^{rank) — c, given as an array rank[0..d — 1] of base-6 digits 

1 forward true / / indicates whether we go through 6rgc''^* forwards or backwards 

2 for i ^ to d — 1 do 



if forward = true then rank[i] ~ c[i] else rank[i] — {b — 1 — c[i]) 
if c[i] is odd then forward -s— -^forward 
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Algorithm 2.4: Comparison operator for a 6-regular curve witli c(r) — 6RGC'^(r) 

Input: Points p — {p[0], ...,p[d — 1]) and q — {q[0], q[d — 1]), witli coordinates in [0, 1), given by 

the digits of the fractional part in hase-b representation 
Output: true if f{p) < f{q), false otherwise 

1 for i^Otod— 1 do reflected[i] ^ false // maintains which columns of T have minus 

2 for i <— to d — 1 do aa:is[i] ^ i // maintains which entry in each row of T is non-zero 

3 initialize an array aHAxis[0..d — 1] // for a new value of axis under construction 

4 repeat 



5 
6 
7 

8 
9 

10 
11 
12 



13 
14 
15 
16 



/ / indicates whether we go through ' forwards or backwards 

qDigit -s— ExTRACTFiRSTDiGiT((j[a2;zs[z]]) 



forward -s— true 
for I ^ to d — 1 do 

pDigit <— EXTRACTFiRSTDiGiT(p[aa;js[i 
if reflected[axis[i]] then 
1^ pDigit <— (b — 1 — pDigit); qDigit -S— (6 — 1 — qDigit) 

if pDigit < qDigit then return forward else if pDigit > qDigit then return -^forward 
if forward then rank[i] ^ pDigit else rank[i] -s— (6 — 1 — pDigit) 
if pDigit is odd then forward -s— -^forward 

for i ^ to d — 1 do 



j -S— ai(rank) 
altAxis[i] 4— axis[j] 

if mj{rank) ~ 1 then refl,ected[altAxis[i]] <~ 

17 swap axis and alt Axis 

18 until aZ? remaining digit strings p[0], — 1] and q[0], 

19 return false 



'.fleeted [ alt A xis [i] ] 



,q[d~ 1] are empty 



1 1 p and q are equal 



3 Proof ingredients 

In this section we establish notation and basic results that are needed for the proofs in the following 
sections. Readers who wish to learn about the final results and pseudocode for implementations only 



without studying the proofs, may skip this section. More precisely, in Section 3.1 we set up notation 
common to most proofs; in Section [3. 2| we set up basic terminology and observations which will allow us 
to reason about interdimensional consistency; in Section |3.3| we do the same for the concept of neutral 
orientation. In Section [X4| we establish basic properties of reflected Gray codes that we will need in the 
interdimensional-consistency proofs; this last section can also be skipped, to be read, if desired, when it 
is referred to later on. 

3.1 Notation for editing numbers, vectors and sequences 

The take-out operator If a is a vector or a number, we will use l^aXo denote the vector or number 
that results from removing the element or digit and renumbering the remaining elements consecutively, 
starting from zero. Thus, if 6 = ^ we have bj = Oj for j < z, and bj = Cj+i for f > i. When the 
take-out operator is applied to a permutation a, the remaining elements are numbered consecutively. 
Thus, we have: 

^ Oj = Oj if < j < i and aj < a.i ; 

^ Oj = flj — 1 if < j < z and Oj > a.i] 

X Oj = flj+i \i i < i < d — 1 and flj+i < a^; 

X Oj = ftj+i — 1 ii i < i < d — 1 and aj^i > a,. 

For example, if a = 30142, then "^a — 2041. The take-out operator can also be applied to sets and 
sequences: if ^4 is a set or sequence, then l^Ais the set or sequence that results from applying ^ to every 
member of the set or sequence. For example, if A is the sequence of ternary numbers (010, 000, 210, 222) , 
then lA\s the sequence (01, 00, 21, 22). 
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The reduce operator If ^ is a set or sequence, then we use '' A to denote the set or sequence that 
results from only selecting those elements a € A such that — j, and applying * to those elements. 
For example, if A is the sequence of ternary numbers (010,000,210,222), then '^^^ A is the sequence 
(01,00,21). 

The reflection operator If a is a base-5 number, then we use q a to denote the number a that results 
from replacing digit by its complement with respect to 5 — 1. For example, if a is the ternary number 
20, then = 22, and if a is the binary number 1001, then qU = 1011. Like the take-out operator, the 
reflection operator can also be applied to sets and sequences. 

The insertion operator If a is a vector or a number, we will use ^^-^ a to denote the vector or number 
that results from inserting an element or digit with value j, and re- indexing the original elements or 
digits Ui, ... as a^+i, .... For example, if a is the ternary number 120, then a = 1210. Like the take-out 
operator, the insertion operator can also be applied to sets and sequences. 



3.2 Visible orders of space-filling curves 

Recall from Definition [2] that an (infinite) set F of space-filling curves is interdimensionally consistent, 
if F contains a unique d-dimensional space- filling curve for any integer d > 1, and fj Cz F is extradi- 
mensional to fi £ F whenever j > i. Observe that, to show that a given set F of space-filling curves is 
interdimensionally consistent, we only need to show that each fj is extradimensional to fj-i. Then it 
follows by induction that fj is also extradimensional to all fi with j > i. 

To make the proofs in the rest of this paper easier to read, we introduce the following terminology. 
The front face F^ of a d-dimensional unit hypercube is the {d — 1) -dimensional face that consists of the 
points {{qo, ■■■,Qd) G [0,1]^^ | Qi = 0}, and the back face F^ of a d- dimensional unit hypercube is the 
{d — 1) -dimensional face that consists of the points {{qo, qd-i) € [0, 1]'' | qi = 1}. With a face of a 
space-filling curve g we mean a face of the unit hypercube filled by g. 

The visible order of a valid ordering y on a face -F* of a d-dimensional space-filling curve g, is the 
function / : [0, l]''^^ [0, 1] that is defined by the following two properties: (i) for any pair of points 
a, 6 e F^ we have: 

7( » < 7( i 6) ^ g{a) < m, 

and (ii) the measure of the set {p \ f{p) G [a, b]} is equal to 6 — a. The second property only serves to 
make / well-defined. The first property is what matters: this property defines the ordering of the points 
^ F^ induced by /. If any valid ordering / of a space-filling curve / can be defined as the visible order 
on F^ of some valid ordering g of g, then we say that g shows f on face _F*. 
Using the above terminology, we get the following: 

Lemma 1 A set F of b-regular space-filling curves, containing a d-dimensional space-filling curve fd for 
any d > 1, is interdimensionally consistent if and only if each curve fd (for d > 1) shows fd-i on each 
front face. 



3.3 Orientations 



Recall from Section 2.3 that for each of the b'^ regions S{r) in a d-dimensional 6-regular mono-curve, a 
permutation a(r) defines the rotation component of the transformation that maps the curve through the 
unit hypercube to the curve in region S{r). When we expand the recursion, we find subregions within 
regions. We can identify each such region by a sequence R in the following way. Let prefix(i?) be R 
without its last element, let last(i?) be the last element of R, and let () be the empty sequence. Then 
S{Q) is the d-dimensional unit hypercube, and S{R) is region S'(last(i?)) within region S'(prefix(i?)). 
The rotation component a{R) of the transformation that maps the curve through the unit hypercube 
to the curve within S{R), is now the permutation a(last(i?)) applied to the permutation a(prefix(i?)), 
where a(()) is the identity permutation. 

We say that a 6-regular mono-curve has neutral orientation, if, in the limit for increasing recursion 
depth i, all d\ permutations are equally frequent among the permutations a{R) of the regions identified 
by sequences R of length £. In other words, if i? is a long random sequence of d-digit base-& numbers, 
then the probability that a{R) is any particular permutation should be 1/dl. 
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Lemma 2 A d- dimensional b-regular mono-curve has neutral orientation if and only if there is a num- 
ber i, such that each of the d\ possible permutations of d numbers can be constructed as the composition 
of £ permutations from the set {a(0), a((6 — 

Proof: Let all possible permutations be numbered from 1 to d\, where permutation 1 is the identity 
permutation. Let P be the dl x dl matrix defined by Pij = zjb'^ if exactly z out of b'^ regions Sir) 
have the permutation that, applied after permutation i, results in permutation j. Note that for each 
permutation a(r) used in z regions S(r\ we put an entry with value z/fc'* in each row i (namely at 
Py , where j is the permutation that results from applying a(r) to permutation i) and in each column j 
(namely at _Py, where i is the permutation that results from applying a(r) to permutation j). Thus each 
row sums up to one and each column sums up to one. Let u be the dl-vector with ux = \ and = for 
ah i > 1. 

Now, if we choose a random region identified by a sequence R of length the probability that it 
has permutation i is given by element i of the vector P^u. As I goes to infinity, P^u converges if and 
only if P is a regular square matrix, that is, if and only if there is a fc such that all elements of P^ are 
strictly positive. Note that all entries in the first column of P^ are strictly positive if and only if each 
permutation can be constructed from the identity permutation as the composition of fc permutations 
from the set {a(0), ...,a((6 — in this case each permutation can actually be constructed from any 

other permutation and therefore all entries in all columns of P^ are strictly positive. 

What is left to prove is that if there is a fc such that all entries of P^ are strictly positive, then 
limf_j.oo -P^w converges to a vector of which all elements have the value l/d!. Note that the vector to 
which lim£_j.oo P^u converges is the unique vector v that satisfies Pv — v. Since all rows of P sum up to 
one, the vector of which all elements are l/d! is a solution to this equation. The lemma follows. □ 



3.4 Properties of reflected Gray codes 



Recall from Section 2.5 that the reflected Gray code 5RGC'' is defined as follows: 5rgc^ is the se- 
quence (0,1,. ..,6 — 1), and &RGC''"'"^, for d > 1, consists of the concatenation of the sequences a ® 
revifodd(a, 6rgc'*) for a = 0, 1, — 1 in order. As noted before, two consecutive elements 6RGC'^(f) 
and 6RGC''(r) differ in only one digit, and the difference in value of this digit is one. In fact, we observe 
that the difference is in the most significant digit in which the base-6 number r differs from f: 

Lemma 3 Let C be a reflected Gray code with base b, indexed by base-b numbers. If Ci{f) ^ Ci{r), then 
i has the smallest value such that fi ^ ri. 

We will now explain some observations about reflected Gray codes with base 2 and 3 that are relevant 
to Sections |4] and |5] of this paper. 

Ternary reflected Gray codes {b — 3) Recall that 6RGC''(r) can be obtained from r by going through 
the digits of r in order from tq to r^-i: whenever a digit is encountered such that = 1, we 'reflect' 
all following digits, that is, we replace by 6 — 1 — rj for all j > i. Note that if 6 = 3, this operation 
leaves the digits 1 unaffected. Thus, 3RGcf(r) = 1 if and only if — 1, and 3RGcf(r) e {0,2} if and 
only ff r, G {0, 2}. 

Lemma 4 (i) If ri is even, then ^ 3RGC''(r) = 3rgc''^^( * r). (ii) If k is even, ^^'^ 3RGC'' = Srgc''^^ 

Proof: Under the above mentioned step-by-step transformation that transforms r into 3RGC'^(r), even 
digits remain even and do not affect other digits. Therefore we have: if G {0, 2} (and therefore, 
3RGcf (r) e {0, 2}), then * 3RGC''(r) = 3rgc''"^( i r). This establishes part (i) of the lemma. 

Furthermore, note that each prefix of 3RGC''(r) only depends on the corresponding prefix of r. There- 
fore, selecting all pairs (r, 3RGC'^(r)) such that 3RGcf(r) = fc, for certain values i € {Q,...,d— 1} and 
fc e {0,2}, results in selecting 6* bundles with consecutive values of r, such that the values of r in 
each bundle are exactly all values that have a common prefix of i + 1 digits, and different bundles 
differ in at least one of the first i digits of r. Thus, replacing each selected pair (r, 3RGC'^(r)) by 
( ^ r, ^ 3RGC''(r)) = ( ^ r, 3rgc'^-1( ^ r)) results in a sequence ordered by ^ r, and thus, the second el- 
ements of the pairs constitute exactly the sequence Srgc''"^ in order. This establishes part (ii) of the 
lemma. □ 
Observe that 3RGC'*((2)d) = (2)^. 
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Binary reflected Gray codes {b = 2) Recall again that 6RGC'^(r) can be obtained from r by going 
through the digits of r in order from vq to r^-i'. whenever a digit is encountered such that = 1, we 
'reflect' all following digits, that is, we replace rj by 6 — 1 — rj for all j > i. For ease of notation, define 
r_i = 0. If = 2, one can verify easily by induction that the end result is the following: 

Lemma 5 2RGcf(r) = jr^ — ri^i\. 

This leads to the following lemma: 

Lemma 6 (i) If Vi = ri_i or i ^ d - 1, then ^ 2RGC''(r) = 2rgc'^^^( ^ r). If i < d - 1 and r; 7^ r^^i, 
then X 2RGC'^(r) = ^2rgc''^^( * r). (ii) For < i < d we have '='^2rgc'^ = 2rgc''^\- forO<i<d-l 
we have 2rgc'^ = i,2RGC'*"^ 

Proof: From Lemma [5] we get that the removal of a single digit does not affect the remaining digits in 
any of the following cases: (a) = r^-i (and thus, 2RGcf (r) = 0); and (b) i = d — 1 (because there are 
no further digits). Thus, in these cases we would also have ^ 2RGC'*(r) = 2rgc''~^( ^ '')• On the other 
hand, if we remove a digit such that i < d — 1 and 7^ r^-i (and thus, 2RGcf (r) 1), the next digit 
of the Gray code changes. Thus we get ^ 2RGC'^(r) = |,2RGC'^~n x f)- 

Following the same reasoning as above in the proof of Lemma Wl selecting all pairs (r, 2RGC'*(r)) such 
that 2RGcf (r) = (for any i) or such that 2RGC^_;^(r) = 1 and replacing each selected pair (r, 2RGC'*(r)) 
by ( X ^, X 2RGC'*(r)) = {^r, 2rgc^~^( * r)) results in the sequence 2rgc''~-^ in order. On the other hand, 
if we remove a digit such that 2RGcf (r) = 1, the next digit of the Gray code changes. Thus, in this case 
we get: * 2RGC'^(r) = q 2rgc''~^( * r); selecting all pairs (r, 2RGC'^(r)) such that 2RGcf (r) = 1 (for any 
i < d - I) and replacing each selected pair (r, 2RGC'^(r')) by ( ^ r, ^ 2RGC''(r)) = ( * r, |, 2rgc'^~^( ^ r)) 
results in the sequence q 2rgc'^~^ in order. □ 
Observe that 2RGC''((l)d) = l(0)d-i. 



4 Extradimensional 3-regular curves 



In this section we study 3-regular vertex-continuous mono-curves, as defined in Section 2.2 We will 
call such curves mono-Wunderlich curves, after Wunderlich |17| . Besides Peano's curve, there are many 
other two-dimensional mono-Wunderlich curves: Figure [7] shows a selection of order-preserving mono- 
Wunderlich curves. Four of these curves we will generalize to higher dimensions: Peano's curve, the 
coil curve [8^ (Luxburg's variation 1 [13]), Luxburg's variation 2 [H [T3] — which we will henceforth call 
half-coil — and the Meurthe curve [5]. 

4.1 Interdimensionally consistent mono-Wunderlich curves 

Each curve described in this section is defined by a single recursive rule. Each curve traverses the S*^ 



subcubes of the unit hypcrcube in the order of the ternary reflected Gray code explained in Section |2.5[ 
that is, for each region S{r) we have: 

c(r) = 3RGCd(r). 

In recursion, this results in each subcube being traversed from one corner to a corner that is opposite 
of it in all dimensions. To make sure that these curve fragments connect up properly, we define the 
reflections m(r) as follows: 

mi{r) = (digitsum(r) — r^) mod 2. 

With these definitions of c{r) and m(r) wc get a proper mono-Wunderlich curve, regardless of our choice 
of the permutations a(r) (Lemma [t] in Section 4.3 proves this). As we will prove in Section 4.4 the 



following choices for a(r) result in interdimensionally consistent sets of mono-Wunderlich curves: 

generalized Peano curves: a(r) is simply the identity permutation {ai{r) = i), regardless of r; 

generalized coil curves: a{r) is the reversal permutation {ai{r) = d — 1 — i), regardless of r; 

generalized half-coil curves: a{r) is the identity permutation if r is odd, and the reversal permuta- 
tion if r is even; 
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Figure 7: Some examples of order-preserving mono-Wunderlich curves. On top we see four curves that 
fit Wunderlich's Serpentine framework, from left to right: Peano, coil, half-coil, and Meurthe. Below, 
from left to right: the 51- or mirrored R-curve (Wunderlich's Meander), and three unnamed curves. 



Table 1: Three-dimensional mono-Wunderlich curves from interdimensionally consistent sets 

rank loc. 



rank 


loc. 


permutation 


inv. 


rcfl. 






Pca/coil/icl/Mt 


Mt 




000 


000 


012/210/210/210 


210 


000 


001 


001 


012/210/012/210 


210 


110 


002 


002 


012/210/210/102 


102 


000 


010 


012 


012/210/012/210 


210 


101 


Oil 


Oil 


012/210/210/210 


210 


oil 


012 


010 


012/210/012/102 


102 


101 


020 


020 


012/210/210/120 


201 


000 


021 


021 


012/210/012/120 


201 


110 


022 


022 


012/210/210/012 


012 


000 


100 


122 


012/210/012/210 


210 


oil 


101 


121 


012/210/210/210 


210 


101 


102 


120 


012/210/012/102 


102 


oil 


110 


110 


012/210/210/210 


210 


110 


111 


111 


012/210/012/210 


210 


000 



exit 



1, 1, 1 
0,0,2 
1,1,3 

2,2 

1, 1 

2,0 

3 

2 

3 

2 

3 



2,0 

1, 1 
2,2 



112 
120 
121 
122 
200 
201 
202 
210 
211 
212 
220 
221 
222 



permutation 
Pea/coil/_icl/Mt 



inv. refi. 
Mt 



exit 



112 
102 
101 
100 
200 
201 
202 
212 
211 
210 
220 
221 
222 



012/210/210/102 
012/210/012/120 
012/210/210/120 
012/210/012/012 
012/210/210/210 
012/210/012/210 
012/210/210/102 
012/210/012/210 
012/210/210/210 
012/210/012/102 
012/210/210/120 
012/210/012/120 
012/210/210/012 



102 110 
201 Oil 
201 101 
012 Oil 
210 000 
210 110 
102 000 
210 101 
210 Oil 
102 101 
201 000 
201 110 
012 000 



1,1,3 
2,0,2 

1,1 

0,0 

1, 1 

0, 

1 

2 

1 

2 

3 

2 

3 



• generalized Meurthe curves: when e {0, 1}, then ai{r) is the number of zeros and ones among 
Ti+i, r(i_i; when = 2, then a,;(r) is d — 1 minus the number of twos in ri^i, ...,rd-i. Alter- 
natively, we can define a(r) through its inverse a(r): the inverse is constructed from the identity 
permutation by moving all indices i such that g {0, 1} to the front and reversing their order. 

Table 14.11 shows the definitions of the three-dimensional versions of each of these curves. The inverse 
permutations are shown explicitly only for the Meurthe curve; for the other curves, each permutation 
is its own inverse. For illustration, Figure |8] shows the graphical definition of the three-dimensional 
Meurthe curve. Of the curves described above, in three and more dimensions only the Meurthe curves 
have neutral orientation — we prove this in Section |4.5| The main results of Section |4] can therefore be 
summarized as follows: 

Theorem 1 The Peano, coil, and half-coil curves constitute interdimensionally consistent sets of order- 
preserving 3-regular vertex- continuous mono-curves. The Meurthe curves constitute an interdimension- 
ally consistent set of order-preserving 3-regular vertex- continuous mono-curves with neutral orientation. 
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Figure 8: The three-dimensional Meurthe curve. 



Inset 2 How the interdimensionaUy consistent niono-Wunderhch curves were found. 

It was not straightforward to guess the permutations for the interdimensionaUy consistent generalizations of the mono- 
Wunderlich curves from the two-dimensional curves, not even for the coil and half-coil curves. In two dimensions there 
are only two numbers to permute and one cannot tell the difference between, for example, a reversal, a left rotation, and 
a right rotation, or the difference between the identity permutation and the permutation that reverses the whole sequence 
except the first element. Even in three dimensions, one cannot tell the difference between a complete reversal and a simple 
swap of the first and the last element. 

To get started, I worked from the assumption that I could require more than what is strictly necessary: I tried to ensure 
that the three-dimensional curves show the two-dimensional curves on all faces instead of only on the front faces. Thus 
it was possible to deduce what the permutations in the eight corner regions had to be, and to some extent, also what the 
permutations in other regions had to be. Having a hypothesis for the necessary permutations in three dimensions, and 
with the experience gained with the harmonious Hilbert curves that are explained in Section [s] it became possible to see 
the patterns emerge. 

Those who are familiar with our previous work in two dimensions may miss two generalizations, namely interdimensionaUy 
consistent generalizations of Wunderlich's Serpentine 011010110 and R-order (Wunderlich's Meander) These curves 
seem to be considerably more difficult to generalize in an interdimensionaUy consistent way than the Peano, coil, half-coil 
and Meurthe curves. The above approach, unfortunately, results in inconsistent requirements on the permutations for the 
three-dimensional generalizations of Serpentine 011010110 and R-order. One would therefore have to aim for a solution 
in which some of the back faces do not show the two-dimensional curve, or, for example, a mirrored version of it. 



4.2 Implementation of comparison operators 

In this section we discuss how to implement comparison operators for each of the four sets of curves 
presented above. We will first explain how to do this for points in the unit hypercube, then explain how 
to extend the operators to points with arbitrary non-negative coordinates, and finally how to extend 
the operators to points with arbitrary real coordinates. The main purpose of this section is to clarify 
how such comparison operators can be implemented, and to enable us to discover more properties of 
the curves by analysing the algorithms. For real implementations, one would have to implement the 
ExtractFirstDigit function, one would have to implement an efficient check to determine if no more 
digits remain to be processed, and one may consider further optimizations, such as look-up tables for 
permutations. 

Modifications for all curves For points in the unit hypercube we can simply fill in the details of 
Algorithm [2^ 

To start with, we substitute "2" for "6 - 1" . 

Second, we implement the reflections. The condition mj(rank) — 1 translates to: rank[j] is odd and 
rank has even digit sum, or rank[j] is even and rank has odd digit sum. Therefore we implement the 
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reflections correctly if we replace the condition on Line 16 by "rank has odd digit sum", and replace 
Line [12] by: 

if pDigit is odd then 
|_ forward <— ^forward; reflected[axis[i\] i— -^reflected[axis[{\] 



Even better, we move Line[5]to just before Line[4] and remove Line 16 altogether. The only effect of this 
is that in some iterations of the main loop the values of forward and all entries of reflected are 'wrong', 
but since all of them are wrong, the 'error' made on Line |9] will be undone on Line 10 or 11 (whichever 
is apphcable). 

A last simplification we make, is that we will work with reflected versions of pDigit and qDigit 
whenever forward is false. This is implemented by changing the condition in Linejsjto {reflected[axis[i]] = 
forward)] changing the return values on Line 10 to true and false, respectively; and changing Line 11 
to simply rank[i] pDigit. This may change the value of pDigit as used on Line |12[ but it does not 
change the parity of pDigit, which is all that matters. 

Thus we get Algorithm |4.1| 



Algorithm 4.1: Framework for implementation of Peano, coil, half-coil and Mcurthc curves 



1 for i to d — 1 do reflected[i] ^ 

2 for i ^ to d — 1 do axis[i] ^ i 

3 initialize an array altAxis[0..d — 1] 

4 forward true 

5 repeat 



false 



/ / maintains which entry in each row of T is non-zero 
// for a new value of axis under construction 



6 
7 

8 
9 

10 
11 
12 
13 



qDigit ^ ExTRACTFlRSTDlGlT((j[aa;is[z]]) 



for z -s— to d — 1 do 

pDigit -It- ExTRACTFiRSTDiGiT(p[aa;js[i] 
if reflected[axis[i]] = forward then 
1^ pDigit (2 — pDigit); qDigit -S— (2 — qDigit) 

if pDigit < qDigit then return true else if pDigit > qDigit then return false 

rank[i] -s— pDigit 
if pDigit is odd then 
1^ forward -S— -^forward; reflected[axis[i]] -^reflected[axis[i]] 



14 for z <— to c? — 1 do altAxis[i] ^ axis[ai{rank)] 

15 swap axis and alt Axis 

16 until all remaining digit strings p[0], — 1] and q[0], . 

17 return false 



,q[d~ 1] are empty 



1 1 p and q are equal 



Peano curves Since ai{rank) = i, regardless of the value of rank, we can simply omit Lines 11 



14 and 15 from Algorithm 4.1 and replace each occurrence of "aa;is[i]" by "i" 



Coil curves Since ai{rank) = d—l — i, regardless of the value of rank, we can implement the updating 
of axis as follows. Omit Lines [TT] and [M] and replace Line [3] by: 

for i to c? — 1 do altAxis[i] (li — 1 — i) 

Half-coil curves In half-coil curves, the permutation is reversed if and only if rank is even, that is, if 
and only if forward remains unchanged. We implement this as follows. Omit Lines [TT] and |14| Before 
Line [6] add: 

direction 4— forward 

Replace Line [TS] by: 

if forward = direction then swap axis and altAxis 
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Meurthe curves Replace Line 14 by: 



indexForNextTwo ^ d; for i <— to d — 1 do if rank[i] = 2 then decrement indexForNextTwo 
indexForNextZeroOrOne -s— indexForNextTwo — 1 
for i to d — 1 do 
if rank[i] — 2 then 
|_ altAxis[indexForNextTwo] aa;is[i]; increment indexForNextTwo 

else rank[i] is or 1 

|_ altAxis\indexForNextZeroOrOne\ -s— aa;is[«]; decrement indexForNextZeroOrOne 



Of course the first for loop can be avoided by integrating it into the loop on Line [6ljT3| of Algorithm |4.1[ 
When implementing this curve, make sure that the algorithm applies a (like the above piece of code 
does), not a. A mistake is easily made and can be hard to detect: as can be seen in Table 
three-dimensional Meurthe curve a(r) equals a(r) for most, but not all, values of r. 



4.1 



in the 



Points with non-negative coordinates We can extend our implementation so that it can also 
compare points outside the unit hypercube. To this end, we should extend the space-filling curve outside 
the unit hypercube. This can be done by considering the unit hypercube as the first subregion of a larger 
cube. The easiest is to consider it as the first subregion of the first subregion, because for all four sets of 
curves, the first subregion of the first subregion is neither reflected nor rotated. Thus, the outcome of a 
comparison between two points that lie in the unit hypercube already, remains the same. Applying this 
idea recursively, we get the following algorithm: 

while any coordinate of p or q is at least 1 do divide all coordinates oi p and g by 9 
run Algorithm |4.1| on p and q 



Points with negative coordinates The above implementation does not work for points with negative 
coordinates. We could try to solve this by considering the unit hypercube to be the subregion in the 
centre of a large cube. However, in that case the entrance point of the larger cube does not lie in the 
origin anymore, and we loose the property that the curve traverses points in an axis-parallel hyperplane 
through the origin in the same order as the lower-dimensional curve. We can solve this as follows. 

Note that, for the coordinates common to all points of an axis-parallel hyperplane through the origin, 
the implementations of the Peano, coil and half-coil curves effectively do not distinguish between the 
digits zero and two: in Line[9j zeros may become twos and vice versa, but as these digits remain even and 
remain the same between p and q, this does not affect the decisions taken by the algorithm, and they do 
not affect forward and reflected. This is also the case for the implementation of the Meurthe curves: a 
digit zero or two for coordinate i only affects where this coordinate ends up between the other coordinates 
in the permutation on Line |14[ but it does not affect the permutation of the other coordinates relative 
to each other. This is all that matters, since coordinate i, being shared by p and g, is not going to decide 
which of p or g comes first in the order in any case. Therefore, all four sets of mono-Wunderlich curves 
presented in this section remain consistent even if the origin of the coordinate system is translated into 
the hypercube by any vector whose coordinates, in ternary representation, only contain zeros and twos. 

Therefore, to allow for negative coordinates, we could consider, for k — 0,1,2,..., each hypercube 
Uk = [\ — jSl'"', \ + |8l'^)'' to be the last subregion of the first subregion of the last subregion of the first 
subregion of the hypercube Uk+i = [\ - 3:8l''+^, \ + |8l''+^)'' (note that the last of the first of the last 
of the first subregion is never reflected nor rotated, in all four sets of curves). Using ternary numbers, we 
write the coordinates of Uk as [jj — jjlOOOO'^, jj -I- 10 000'^)'*. We will continue to use ternary numbers 
in the rest of this paragraph until further notice. Define So{x) = x and let Si{x) = (x-f 202)/10 000 be the 
function that maps the coordinate range [jj — ytIOOOO'^, jj + yjIOOOO'^] of Uk to the coordinate range 

- j^lOOOO'=-i, ^ + iflOOOO'^-i] of Uk-i- Let Sk{x), for fc > 1, be given by Sk{x) = si{sk-iix)). 
We can now define a local coordinate system for each hypercube Uk, were Sk{x) gives the value in the 
local coordinate system that corresponds to x in the global coordinate system. The origin of the unit 
hypercube Uq now lies at position (sfc(O), 5^(0), ...) in the local coordinate system of hypercube Uk- Note 
that s/c(0) is the number 0.(0202)^. Thus, the origin lies at coordinates that only contain zeros and 
twos in the local coordinate system of each hypercube Uk , and hence, the space-filling curve filling Uk is 
extradimensional to the lower-dimensional space-filling curves with respect to this origin. 
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The following algorithm (with numbers in decimal format again) zooms out to the smallest hypercube 
Uk that contains the points to be compared, while converting the coordinates to the local coordinate 
system of Ut- 

while any coordinate of p or q is less than zero or more than one do 

add 20 to all coordinates of p and q // adding the ternary number 202 

divide all coordinates of p and 5 by 81 // dividing by the ternary number 10 000 

run Algorithm |4.1| on p and q 

Note that this solution produces different results than the solution for points with non-negative coordi- 
nates that was presented before. 

4.3 Proof of vertex-continuity 

Lemma 7 The generalizations of the Peano, coil, half-coil and Meurthe curves defined above, indeed 
constitute mono-Wunderlich curves. 

Proof; All we have to prove is that the curves are vertex-continuous. 

Recall that the curves start in the region with location (0)^, and end in the region with location (2)^. 
Thus, for the first and the last subregion r, we have mi{r) = for all < i < d. As a result, in the 
recursion within S'(O), the curve will always start in the region closest to the origin, and in the recursion 
within S{{2)d)^ the curve will always end in the region farthest from the origin. Thus the entrance gate 
of the curve is at the origin, and the exit gate is at the corner of the c?-dimensional hypercube that is 
farthest from the origin. 

Thus, for the entrance gate of a subregion S{r) to be connected to the exit gate of the previous 
subregion S'(f), region S{r) should be mirrored in all dimensions with respect to 5'(f), except in the one 
dimension in which the locations of S{r) and S{f) differ. More precisely, (i) if Ci{f) < Ci{r), we should 
have mi{f) — mi{r) — 0; (ii) if Ci(f) = Ci(r), we should have mi{f) + mi{r) = 1, and (iii) if Ci(f) > Ci(r), 
we should have mi{f) ~ mi{r) — 1. 

First consider case (i): Ci{f) < Cj(r). Since Cj(r) ^ Ci(r), we have that i has the smallest value such 
that f-i ^ Ti, and therefore, ri = fi + 1. We have fo...ri_i = r^.-Xi-i. Moreover, since Ci(f) < Ci(r), we 
have, by definition, 3RGCi(f) < 3RGCi(r), and the number of reflections caused by the digits tq, ...,ri_i, 
that is, the number of odd digits among 3RGCo(r), 3RGCi_i(r), must be even. Since 3RGCj(r) is odd 
if and only if rj is odd, for any j e {0, — 1}, this implies that the number of odd digits among 
ro,...,ri_i must be even. Therefore, X]}=o ~ X]j"=o = (mod 2). By definition, ?' = f + 1, and 
therefore, since r and f are ternary numbers, we have ri+i...fd_i = {2)d-i-i and ri+i...rd^i = (0)d_i^i, 
so = ''j = (mod 2). Recall that mi{r) is defined as (digitsum(r) — r^) mod 2, so we 

have m,(f) = (E}=o + J2pl+i ^j) mod 2 = 0= (E}=o + Ej^^+i ^j) mod 2 = m,{r), QED. 

In case (ii), Ci{f) = Cj(r), we know that i is not the smallest value j such that rj 7^ r^. We have 
either i < j and ri = fi, or i > j and = 2 and = 0; in both cases = fi (mod 2). Note 
that for ternary numbers, r = f + 1 implies digitsum(r) = digitsum(r) -I- 1 (mod 2). Thus we have 
mi{r) = digitsum(r) — r.; = digitsum(r) + 1 — ri = digitsum(f) + 1 — f; = nii{f) + 1 (mod 2), which 
implies mi{r) = {rriiif) + 1) mod 2, QED. 

Third, consider case (iii): Ci(f) > Ci{r). The proof is analogous to case (i), except that now, since 
Ci{f) > Ci{r), the number of odd digits among ro...ri_i must be odd, and thus we get mi{f) — (X]j=o 
YfjZ^i+i fj) mod 2 = 1 = (Ej=o ^3 + TfjZl+i rj) mod 2 = m,(r), QED. □ 

4.4 Proof of interdimensional consistency 

In this section we will prove that the above mentioned generalizations of the Peano, coil, half-coil and 
Meurthe curves are interdimensionally consistent. Arguments to this effect were also given in the para- 
graph of negative coordinates in Section |4.2[ However, those arguments depend on the specific tie- 
breaking rule for points that lie on the boundary between regions, that is, the choice of /, that is 
implemented in the pseudocode in this paper. Moreover, we need to develop a different proof technique, 
because in Section |5] a look at the pseudocode will not reveal the secrets of the curves as easily. There- 
fore, in the present section, we will prove interdimensional consistency directly from the definitions of 



20 



the curves, rather than from the implementation in pseudocode. I will first describe the part of the proof 
that is common to all four generalizations, and after that we will discuss the details of each separate 
generalization. 

Consider any of the four sets of space-filling curves defined above, containing a d-dimensional space- 
filling curve fd for any d > 1. As observed in Lemma [T| to prove that the set is interdimensionally 
consistent, it is sufficient to show that each curve fd (for d > 1) shows fd-i on each front face. In fact, 
we will prove something stronger, namely that each curve fd shows fd-i on each face. 

Below, we write / for fd, and by c, a, m, o, M, and r we denote the location, permutation, reflection, 
translation, transformation matrix and transformation functions that define /. We write /' for fd-i, and 
by c', a', m', o', M' and r' we denote the location, permutation, reflection, translation, transformation 
matrix and transformation functions that define /'. Our proof that each curve / shows /' on each face 
goes by induction on increasing level of refinement. More precisely, for a level of refinement £, we consider 
the recursive balanced subdivision of the (rf— l)-dimensional unit hypercube U' into 3^*^'*~^-' regions, and 
prove the following induction hypothesis: the visible order of / on any face F' = visits these regions 
in the same order as /'. 

Base case 1=1 As a base case, consider refinement level £ = 1, and a face F' = F^ . Recall that 
/ is defined based on subdivision of a d-dimensional unit hypercube — which we will denote U — into 3*^ 
regions S'(O), S{{2)d). This induces a subdivision of F' into i'^''^ regions F'(0), F'{{2)d-i), each of 
which is a (d — 1) -dimensional face of one of the regions S'(O), S'((2)d). Let F'(0), F'((2)(i_i) be 
indexed in the order in which the c?-dimensional regions that contain them appear in 5'(0), S{{2)d). 
Thus, there is a monotonously increasing function a : {0, {2)d-i] — {0, (2);;} such that F'{r') is a 
face of S{a(r')). Note that for r' e {0, (2)^-1}, the locations of the regions F'ir') within U' are given 
by X c(iT(r')) = X 3RGC''((T(r')). Furthermore, the co-domain of a consists of exactly those values r such 
that Ci{r) ~ 3RGcf (r) = 2k. Thus, the locations of the regions F'{r') within [/', in order of increasing r', 
form the sequence ^^•^'^ 3rgc'^, which, by Lemma |4j is exactly the sequence 3rgc''~^ = c'. Hence, the 
visible order of / on any face F' = F^ visits the regions of refinement level ^ = 1 in the same order as /'. 

Induction We will now prove that if the induction hypothesis holds for refinement level f — 1, then it 
also holds for refinement level £. Again, consider a face F' = F^ , subdivided into 3*^"^ top-level regions 
i^'(0),...,F'((2)d-i), each subdivided into smaller regions. By the analysis of the base case, 

the S''"^ top-level regions i^'(O), F'((2)d_i) of F' are traversed in the correct order. What is left to 
prove, is that for each r' G {0, (2)d-i}, the visible order of the region F'{r') traverses its recursive 
subdivision into 3^^^^-'^'*^^-' smaller regions in the correct order. 

The transformation T(r) that transforms the curve / through U into the section of / through S{a{r')), 
also transforms a face F oiU into F'{r'). Denote this face F by F*{r'), denote ai{a{r')) by j, and observe 
that we have: 

Note that by the inductive assumption F*{r') shows the correct order /' of its subregions of level I — \ 
(which map to subregions of level £ on F')-. all that is left to prove is that the transformation T(cr(r')) 
that operates on S{a{r')) transforms /' on F*{r') correctly into /' on F'ir'), that is, we need to show 

l,T{a{r')){^=^^U')=T'{r'){U'), 

where ^^^^ U' is the operation that inserts a coordinate pj = 2k in the coordinate sequence of each point 
p e U' . Recall that T(r) is the function p — >• ^M{r)p + o{r), so what we have to prove is: 

I (iM(a(r'))(^f '^p) + lo{o{r'))) = \m' {r')p + o' {r'). 

By the definition of s(r') the translations o{a{r')) and o'(r') are such that F*{r') is mapped to F'{r'), so 
we do not need to worry about those. Our challenge is to show that the matrices M{(T{r')) and M'(r') 
match, that is, M'{r') must equal M{a{r')) with row i and column j removed. Given the definition of 
these matrices, we can separate this challenge into two subproblems. First, the reflections should match: 

;,m{air'))=m'{r'), (1) 
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and second, the permutations should match: 

Xa(/))=a'(r'). (2) 

We start with the reflections. To prove is Equation [l] Let r be cr{r'). By definition of r — cr(r'), we 
have c(r) = *^^''c(r'), and hence, 3RGC,f(r) = 2k and ^ 3RGC''(r) = 3RGC'^~^(r'). Recall that 3RGcf(r) 
is even if and only if is even. Therefore, is even, and by Lemma |4] part (i), we have: 

r'-^r=Xr'), (3) 

and: 

r' = r = a{r') (mod 2). (4) 

By the definition of the curves, ni'^{r') = digitsum(r') — (mod 2), and because r' is a ternary number, 
we have digitsum(r') = r' (mod 2). Now we have, for j < i, m'j{r') = digitsum(r') — = r' — r'^ = 
cf{r')~Uj{r') = mj{a{r')) (mod 2), and for j > i, m'j{r') = digitsum(r')— = r'—r'j = a{r')—(jj+i[r') = 
mj+i{cj{r')) (mod 2). Thus m'(/) = * m(CT(r')), QED. 

Now all that is left to prove is that the permutation components are correct, that is, we need to prove 
Equation [2j Below we prove this separately for each of the four sets of curves described above. 

Peano curves Since a'(r') and a{a{r')) maintain the order of all coordinate axes regardless of r' and 
cr(r'), we always have a'(r') — * a{a{r')), QED. 

Coil curves Since a'(r') and a{a{r')) reverse the order of all axes regardless of r' and cr(r'), we always 
have a'{r') = ^ a{a{r')), QED. 

Half-coil curves The permutation a'(r') reverses the order of all axes if r is odd, and maintains the 
order of all axes if r' is even. The permutation a{ij{r')) reverses the order of all axes if cr(r') is odd, and 
maintains the order of all axes if cr(r') is even. By Equation |4] we have that r' is odd if and only if (T(r') 
is odd, and thus, the permutations a'(r') and o(a(r')) match, that is, a'(r') — ^ a{a{r')), QED. 

Meurthe curves We will again use Equation |3] Recall that the respective inverses of a'(r') and a{r) 
(where r — <j{r')) can be constructed from the identity permutation by moving all indices h such that 
I'hi 01' respectively, is or 1 to the front, and reversing the order of these indices. Obviously, inserting 
a digit in r' at position i to obtain r — cr(r'), moving the index i in a(r), and subsequently removing 
index i from a(r) (which is what ^ a('') does), does not affect the order of the remaining indices relative 
to each other. Thus wc have a'(r') ~ ^ '^{'^{''''))^ QED- 

4.5 Orientation of the interdimensionally consistent mono-Wunderlich curves 

By Lemma [2j a d-dimensional mono-Wundcrlich curve has neutral orientation if and only if there is 
a number H., such that each of the d\ possible permutations of d numbers can be constructed as the 
composition oil permutations from a. Trivially, each 1-dimensional curve has neutral orientation. Below 
we analyse whether the d-dimensional curves for d > 2 have neutral orientation. 

Peano curves In the Peano curves, a only contains the identity permutation; therefore no other 
permutation can be constructed and the Peano curves do not have neutral orientation. 

Coil curves In the coil curves, a only contains the reversal permutation. As a result, there is no I 
such that each permutation can be constructed as the composition oi I permutations from a, not even if 
d — 2: if is odd, only the reversal permutation can be constructed, and if i is even, only the identity 
permutation can be constructed. Therefore the coil curves do not have neutral orientation. 

Half-coil curves In the half-coil curves, a only contains the identity and the reversal permutation. As 
a result, there is no £ such that each permutation can be constructed as the composition of i permutations 
from a, unless d — 2. Thus, the two-dimensional half-coil curve has neutral orientation, but the higher- 
dimensional half-coil curves have not. 
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Meurthe curves In the Meurthe curves, a(00(2)£;_2) is the permutation that swaps the first two 
elements, and a((2)d-i0) is the permutation that rotates the whole sequence one step to the right. 
Thus, we can swap any two adjacent elements by composing at most d + 1 permutations: d rotations 
a{{2)a_i0) and one swap a(00(2)d_2), placed between the rotations at the point where the two elements 
to be swapped have been rotated to the first two positions. Since any permutation can be created 
from any other by at most d{d — l)/2 swaps of adjacent elements, we get that any permutation can be 
created from any other by composing at most {d+l)d{d— l)/2 permutations a{00{2)d-2) and a((2)rf_i0). 
However, to satisfy the condition in Lemma[2] there should be a depth £ such that each permutation can 
be created from the identity permutation by composing exactly £ permutations. In the Meurthe curves, 
this can indeed be done for £ = (d + l)d{d — l)/2, since any shorter sequence of permutations can be 
completed by the identity permutation found at a((2)d). If follows that the Meurthe curves have neutral 
orientation. 

Note that the analysis given above is by no means tight: in reality much less than {d + l)d{d — l)/2 
permutations from a may suffice to produce all possible permutations. For example, in four dimensions 
two permutations suffice, rather than [d + l)d{d — l)/2 = 30. 

5 Extradimensional 2-regular curves 

In this section we study 2-regular vertex-continuous mono-curves, as defined in Section [2?2l We will call 
such curves mono-Hilbert curves [5]. In two dimensions, the properties that define a mono-Hilbert curve 
form a minimal sufficient set of properties that defines Hilbert's original curve (see Figure [6|. Therefore 
one can argue, although it is ultimately a matter of taste, that in higher dimensions these properties are 
necessary and sufficient for a curve to be considered a true generalization of a Hilbert curve [5]. Below 
we will discuss two sets of mono-Hilbert curves in arbitrary dimensions: the generalization implemented 
by Moore 1141 (citing Butz [3]), and a new generalization called harmonious Hilbert curves. 

5.1 Standard Hilbert curves 

We will first describe a class of Hilbert curves in higher dimensions, in which the order of the subregions 
and the connection points between them follow a very regular pattern. We call these curves standard 
Hilbert curves (by lack of a better name). 

Each standard Hilbert curve is defined by a single recursive rule. Each curve traverses the 2^* subcubes 
of the unit hypercube in the order of the binary reflected Gray code explained in Section [2. 5[ that is, for 
each region S{r) we have: 

c(r) = 2RGCd(r). 

Standard Hilbert curves differ in their permutation functions a, but oq is fixed, as follows (recall that we 
write f for r — 1, and we write (r) for r + I): 

ad-i{r) =0 if r = (the first subregion) or r = (1)^ (the last subregion); 
ai(r) = if < r < (1)^, r is odd and Ci{r) ^ Ci{f); 
ai{r) — if < r < {l)d, r is even and Ci{r) ^ Ci{f). 

Reflections are defined as follows: 

mi{r) = Cj(r) if r = 0; 

niiir) — \ — Ci{r) if r > and i = d — \] 

mi{r) = Ci{f) if r > and i < d — 1. 

We consider two sets of standard Hilbert curves: 

• Butz-Moore curves, as implemented by Moore [13], based on the work of Butz [3]: here aj{r) — 
(ao{r) +j) mod d. 

• our new harmonious Hilbert curves: when ^ rd-i, then ai{r) is the number of digits among 
ri_|_i, ...,rd^i that differ from r^-i', when rj = r^-i, then ai{r) is d — 1 minus the number of digits 
among tq, ...,ri_i that equal r^-i- Alternatively, we can define a(r) through its inverse a(r): the 
inverse is constructed from the identity permutation by moving all indices i such that ^ r^-i to 
the front and reversing their order, and reversing the order of the remaining indices as well. 



23 



Figure 9: The three-dimensional harmonious Hilbert curve. Left: definition. Right: location of the gates 
between the subregions. 



Table 2: Butz-Moore curves and harmonious Hilbert curves in five dimensions. 



rank 


loc. 


00000 


00000 


00001 


00001 


00010 


00011 


00011 


00010 


00100 


00110 


00101 


00111 


00110 


00101 


00111 


00100 


01000 


01100 


01001 


01101 


01010 


01111 


01011 


OHIO 


01100 


01010 


01101 


01011 


OHIO 


01001 


01111 


01000 



permutat. 
B-M / HH 
12340/43210 
23401/32104 
23401/43201 
34012/21043 
34012/43021 
23401/21403 
23401/43102 
40123/10432 
40123/40321 
23401/24103 
23401/41302 
34012/14032 
34012/41032 
23401/14302 
23401/42103 
01234/04321 



mv. perm. 
B-M / HH 
40123/43210 
34012/32104 
34012/34210 
23401/21043 
23401/24310 
34012/31042 
34012/32410 
12340/10432 
12340/14320 
34012/32041 
34012/31420 
23401/20431 
23401/21430 
34012/30421 
34012/32140 
01234/04321 



refl. 



00000 i 
00000 I 
00000 I 
00011 I 
00011 
00110 
00110 
00101 
00101 
01100 i 
01100 
01111 
01111 
01010 
01010 i 
01001 I 



exit 

0,0,0,0,1 
0,0,0,1,1 
0,0,0,2,1 
0,0,1,2,1 
0,0,2,2,1 
0,0,2,1,1 
0,0,2,0,1 
0,1,2,0,1 
0,2,2,0,1 
0,2,2,1,1 
0,2,2,2,1 
0,2,1,2,1 
0,2,0,2,1 
0,2,0,1,1 
0,2,0,0,1 
1,2,0,0,1 



rank 

10000 
10001 
10010 
10011 
10100 
10101 
10110 
10111 
11000 
11001 
11010 
11011 
11100 
11101 
11110 

11111 



loc. 

11000 
11001 
11011 
11010 
11110 

11111 

11101 
11100 
10100 
10101 
10111 
10110 
10010 
10011 
10001 
10000 



permutat. 
B-M / HH 
01234/04321 
23401/42103 
23401/14302 
34012/41032 
34012/14032 
23401/41302 
23401/24103 
40123/40321 
40123/10432 
23401/43102 
23401/21403 
34012/43021 
34012/21043 
23401/43201 
23401/32104 
12340/43210 



mv. perm. 
B-M / HH 
01234/04321 
34012/32140 
34012/30421 
23401/21430 
23401/20431 
34012/31420 
34012/32041 
12340/14320 
12340/10432 
34012/32410 
34012/31042 
23401/24310 
23401/21043 
34012/34210 
34012/32104 
40123/43210 



refl. 



01001 h 

11000 I 

11000 I 
11011 I 
11011 I 
11110 I 
11110 I 
11101 I 
11101 I 
10100 I 
10100 I 
10111 I 
10111 I 
10010 I 
10010 I 
10001 I 



exit 

2,2,0,0,1 
2,2,0, 1, 1 
2,2,0,2,1 
2,2,1,2,1 
2,2,2,2, 1 
2,2,2, 1, 1 
2,2,2,0, 1 
2, 1,2,0, 1 
2,0,2,0, 1 
2,0,2, 1, 1 
2,0,2,2, 1 
2,0,1,2,1 
2,0,0,2,1 
2,0,0,1,1 
2,0,0,0, 1 
2,0,0,0,0 



Table ISTl shows the definitions of the five-dimensional versions of each of these curves. For illustra- 
tion, Figure [9] shows the graphical definition of the three-dimensional harmonious Hilbert curve. The 
harmonious Hilbert curves are interdimensionally consistent and have neutral orientation (as we will see 
below), the Butz-Moore curves are not interdimensionally consistent (as observed in Section [ij and do 
not have neutral orientation. The main result of Section [5] can be summarized as follows: 

Theorem 2 The harmonious Hilbert curves constitute an interdimensionally consistent set of order- 
preserving 2-regular vertex- continuous mono-curves with neutral orientation. 



5.2 Implementation of comparison operators 

In this section we discuss how to implement comparison operators for each of the two sets of curves 
presented above. We will first explain how to do this for points in the unit hypercube, then explain how 
to extend the operators to points with arbitrary non-negative coordinates, and finally how to extend 
the operators to points with arbitrary real coordinates. The main purpose of this section is to clarify 
how such comparison operators can be implemented, and to enable us to discover more properties of 
the curves by analysing the algorithms. For real implementations, one would have to implement the 
ExtractFirstDigit function, and one may consider optimizations such as look-up tables for permuta- 
tions. You may find that many steps of the algorithm can be implemented and quite easily and efficiently 
by bit operations such as bit shifting and exclusive or. Readers who are not interested in implementations 



may want to skip to Section 5.3 right away. 
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Inset 3 How the harmonious Hilbert curves were found. 

The first non-trivial harmonious Hilbert curve is the three-dimensional curve. This curve was discovered when generating 
all 10,694,807 different three-dimensional mono-Hilbert curves, including curves that are not order-preserving [6]. In fact, it 
turned out that in three dimensions, the harmonious Hilbert curve, showing the two-dimensional Hilbert curve on five faces, 
is the only mono-Hilbert curve that is extradimensional to the two-dimensional Hilbert curve. There is no three-dimensional 
mono-Hilbert curve that shows the two-dimensional Hilbert curve on all six faces (even when allowing reflections on some 
faces). 

Having discovered the three-dimensional curve, the question arose whether similar curves would exist in higher dimensions, 
and if so, how many of them. Unfortunately, the approach to finding interdimensionally consistent generalizations of 
mono-Wunderlich curves, as explained in Inset [2] could not be applied to mono-Hilbert curves: it would already fail in 
three dimensions. Instead, I programmed an exhaustive search algorithm that, given a d-dimensional curve would 
search all possible symmetric {d + l)-dimensional curves that traverse the subregions in binary reflected Gray code order, 
reporting those curves that are extradimensional to /^j. Using effective pruning, unique solutions for four-, flve- and six- 
dimensional curves were quickly found. One may regard this as a somewhat independent confirmation that the proofs given 
in Section |5.4| are likely to be correct. 

It became clear that the challenge was not to run the exhaustive search, but to understand the structure in the solutions. 
Freek van Walderveen was the first to recognize structure in the permutation patterns, and managed to write an efficient 
algorithm that correctly predicted the permutations of the seven-dimensional curve, which was subsequently verified by 
my exhaustive search algorithm. Insights that led to a simpler algorithm began to come when I drew a graph of the 
subregions whose permutations differ by only one inversion. After drawing the graph for the five-dimensional curve on a 
geometric drawing of the five-dimensional hypercube, it became clear that inversions in the permutations of subregions 
S{r) correspond to inversions of adjacent digits with different values in the binary representations of the ranks r. With an 
analysis of the necessary base cases, this led to the simple definition given in Section |5.1| which specifies the permutation 
for any given rank. 



Modifications for all curves For points in the unit hypercube we can simply fill in the details of 



Algorithm To start with, we substitute "1" for "6- 1". 

Second, we implement the reflections. First, consider the reflections of regions with even rank r. Since 
r is even, we either have r = 0, or there is some i such that fi-ifi...fci-i = 0{l)d~i and ri-iri...rd~-i — 
l(0)(j_i, while the first i — 1 digits of r and r are equal. By Lemma [5j i is now the largest i such that 
Ci{r) = 1. Furthermore, by Lemma [sj we have Ci_i(r) Ci_i(r), and thus, by definition, ao{r) = i — 1. 

If r is odd, we either have r = {l)d, or there is some i such that ri-iri...rd-i = 0{l)d-i and 
fi^ifi...fd-i = l(0)d-i, and we find that i is the largest i such that Ci{r) = 1, Ci_i(r) ^ Ci_i(f), and 
ao{r) = i — 1. Finally, observe that r is odd whenever the number of ones in c(r) is odd, since c(0) = 
and, as r increases to (l)c;, the value of r and the number of ones in c(r) both change by one in every 
step. 

This leads to the following simple method to determine the reflections m(r) based on c(r), for < 
r < {l)d. Go through the digits Ci{r) in order of increasing i: if Ci{r) = 0, set mi{r) — 0; if Ci(r) = 1, 
set mi{r) = 1 and set ao{r) = i — 1 (mod d). This procedure computes TOi(r) correctly except when 
Ci{r) Ci(f) OT i < d — 1. Therefore we correct m(r) by changing md-i{r) and, when the number of ones 
in c(r) is even, changing mi{r) for i = ao(r). For r = 0, the exact same procedure also works correctly if 
ao(r) is initialized to c? — 1, as aQ{r) will remain unchanged, mi{r) will be set to zero correctly for all i, 
and in the end, the two corrections (of md-i{r) and of mi{r) for i = ao(r) = d — 1) cancel each other's 
effect. For r = (1)^, we have c(r) ~ l(0)rf_i, and the procedure will also set ao(r) and m,- correctly. 

Therefore we implement the reflections correctly if we modify Algorithm |2.4| as follows. After Line [5) 
add: 

iMinusOne -s— d — 1; altAxis[Q\ ^ axis [iMinus One] 



Replace Line 12 by 



if pDigit — 1 then 

forward -s— -^forward; reflected[axis[i]] ^ -^reflected[axis[{\ 
altAxis[0\ -i— axis[iMinusOne] 

iMinusOne 4— i 
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Remove Line 16 Instead, add the following before Line 17 (outside the for loop): 

reflected[axis[d — 1]] -s— -^reflected[axis[d — 1]] 
if forward then 

// pDigit was one an even number of times 
reflected[altAxis[0]\ -s— -^reflected[altAxis[0]] 

Even better, we move Linejsjto just before Line[4j and 'correct' reflected[altAxis[0]] unconditionally. The 
only effect of this is that we will enter the next iteration of the repeat... until loop with wrong values 
for forward and refiected[axis[0]], but the reader may verify that this does not change the return values 
and rank[i] as computed in the first iteration of the large for loop. In contrast, the value of pDigit will 
change, and thus the values of forward and refiected[axis[0\] will be correct again by the end of the first 



iteration of the for loop. Thus we get Algorithm 5.1 



Algorithm 5.1: Framework for implementation of Butz- Moore and harmonious Hilbert curves 



1 for i to d — 1 do reflected[i] < 

2 for i <~ to d — 1 do aa;i,s[i] ^ i 

3 initialize an array altAxis[0..d — 1] 

4 forward -S— true 

5 repeat 

6 
7 

8 
9 
10 



false 



/ / maintains which entry in each row of T is non-zero 
/ / for a new value of axis under construction 



11 

12 
13 
14 
15 

16 



axis[iMinusOne] 

qDigit 



EXTRACTFlRSTDlGIT((j[a2;zs[i]]) 



iMinusOne ^ d — 1; altAxis[0] 
for z to d — 1 do 

pDigit <— EXTRACTFiRSTDiGiT(p[aa;js[i] 
if refiected[axis[i]] then 
1^ pDigit (1 — pDigit); qDigit -S— (1 — qDigit) 

if pDigit < qDigit then return forward else if pDigit > qDigit then return -^forward 
if forward then rank[i] ^ pDigit else rank[i] (1 — pDigit) 
if pDigit = 1 then 

forward ^ -^forward; refiected[axis[i]\ 4— -^reflected[axis[i]] 
altAxis[0] -i— axis [iMinusOne] 

iMinusOne -f- i 



17 for i -s— 1 to d — 1 do altAxis[i] axis[a 

18 refiected[axis[d — 1]] -^reflected[axis[d - 

19 swap axis and altAxis 

20 until all remaining digit strings p[0], — 1] and q[0\, 

21 return false 



rank)] 

1]]; refiected[altAxis[0]] 



efiected[altAxis[0] 



,q[d— 1] are empty 



1 1 p and q are equal 



Butz-Moore curves The permutations of Butz-Moore curves are always rotations (in the permutation 
sense of the word), therefore they can be maintained by only keeping track of a2:zs[0], without storing 
axis\V..d— 1]. Thus we can simplify Algorithm 5.1 substantially and get Algorithm 5.2 here we maintain 
aa;is[0] in a variable called rotation. We define the value of —1 mod d as d — 1. 



Harmonious Hilbert curves To implement the harmonious Hilbert curves, remove Lines [61 [T5|and[T6l 
and replace Line [17] by Algorithm [53] The first for loop of Algorithm [5 . 3[ can be avoided by integrating 
it into the loop on Line]?- 16 of Algorithm 5.1 



Points with non-negative coordinates To enable our implementation to compare points outside 
the unit hypercube, we consider the unit hypercube as the first subregion of a larger cube. Again, as 
in Section [4?2j for the harmonious Hilbert curves it is easiest to consider the unit hypercube as the first 
subregion of the first subregion of a larger cube, because the first subregion of the first subregion is 
neither rotated nor refiected. Thus we get the following algorithm: 

while any coordinate of p or q is at least 1 do divide all coordinates of p and g by 4 
run Algorithm 5.1 for harmonious Hilbert curves on p and q 
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Algorithm 5.2: Implementation of the Butz-Moore curves 



1 forward 

2 repeat 



true; rotation -s— 0; for 



to — 1 do reflected[i] -s— false 
(i — 1) mod d; newRotation ^ iMinusOne 

EXTRACTFlRSTDlGIT((7[i]) 



i rotation; iMinusOne 
repeat 

pDigit -S— ExTRACTFlRSTDlGlT(p[j]); qDigit 
if reflected[i] then 
1^ pDigit -s— (1 — pDigit); qDigit ^ (1 — qDigit) 

if pDigit < qDigit then return forward; else if pDigit > qDigit then return 
if pDigit — 1 then 
1^ forward -S— -^forward; reflected[i] -S— -^refiected[i] 

iMinusOne ^ i; i ^ (i + 1) mod d 
until i = rotation 

refiected[iMinusOne] -s— -^refiected[iMinusOne] 
refiected[newRotation] -s— -^reflected[newRotation] 
rotation 4— newRotation 

16 until a/Z remaining digit strings p[0], — 1] and ^[0], .. 

17 return false 



■^forward 



newRotation -s— iMinusOne 



,q[d— 1] are empty 



Algorithm 5.3: Algorithm to compute the inverse permutation of region rank of a harmonious 
Hilbert curve, and to apply it to the permutation axis to obtain altAxis. 

1 numberOfOnes 0; for « ^ to d — 1 do numberOfOnes -f- numberOfOnes + rank[i] 

2 if ranA; [d — 1] = then 

3 1^ indexForNext[l] -S— 0; mdea;ForA'^ea;i[0] -S— numberOfOnes 

4 else ranA; [d — 1] = 1 

5 1^ mdexForA'^exip] 0; indexForNext[l] d — numberOfOnes 

6 for i ^ d — I down to do 

7 I altAxis[indexForNext[rank[i]]] 4— a2;is[i]; increment indexForNext[rank[i]] 



Points with negative coordinates One could try to extend the implementation to negative coordi- 
nates in a way similar to what is explained in Section |4j2j Recall that the basic idea is to move the origin 
of the coordinate system for the original, hypercube-filling curves, to somewhere properly inside the unit 
hypercubes. By 'zooming out', one can then obtain curves that fill the full d-dimensional space, rather 
than covering only points with non-negative coordinates. However, the resulting curves would only be 
interdimensionally consistent, if the original hypercube-filling curves retain interdimensional consistency 
with the redefined origin. Since I have not found any location for the origin that meets this requirement, 
I omit the details for extensions to negative coordinates. 

5.3 Proof of vertex-continuity 

In this section we will prove that the Butz-Moore and the harmonious Hilbert curves are indeed curves, 
that is, they are vertex-continuous. We will do so in two steps: first we prove that standard Hilbert 
curves are vertex-continuous, and then we proof that the Butz-Moore and the harmonious Hilbert curves 
are standard Hilbert curves. 

Lemma 8 Standard Hilbert curves constitute mono-Hilbert curves. 

Proof: All we have to prove is that the curves are vertex-continuous. 

We first establish the location of the entrance and exit gates of the curves. Recall that the curves 
start in the region 5(0) with location (0)d, and by definition mi(0) = for all < i < d. As a result, 
in the recursion within iS'(O), the curve will always start in the region closest to the origin. Thus the 
entrance gate of the curve is at the origin. Let r be 2'' — 1 = (l)^;. The curves end in the region S{r) with 
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location c(r) = 1(0)^-1, that is, the region touching the point (1,0, ...,0). Note that c(f) — l(0)£i_2l, 
and therefore, m(r) = l(0)(j-2l- The permutation a(r) is the reversal permutation. The coordinates 
of the last subregion within S{r) are obtained by applying the permutation specified by 0(7") and the 
reflections specified by m(r) to the coordinates c(r): we find that the last subregion visited within S{r) 
is again the region closest to the the point (1,0, ...,0). Thus, in recursion, the curve will always end in 
the region touching the point (1, 0, 0), which is therefore the exit gate of the curve. 

I claim that for < r < 2^*, the connecting gate between subregion S{f) and S{r) is at the point 
x{r) = (xo(r), Xd-i{r)), where Xi{r) = |(ci(r) + Ci{r)) for < i < d—1, and Xd-i{r) = \. For ease of 
notation, we let x(0) = (0, 0, .., 0) and x{2 ) = (1, 0, 0) be the entrance and the exit gate, respectively, 
of the complete curve. (Figure |9] shows an example in three dimensions). To prove that my claim is 
correct, we should establish that the curve within each region S{r) starts in x{r) and ends in x{f). I leave 
it to the reader to verify this for the special cases of the first subregion (r = 0) and the last subregion 
(r = 2"^ — 1). We will now see how to verify that the curve within S{r) starts in x{r) and ends in x{f) 
for < r < 2^* - 1. 

Observe that the location of S{r) differs from the location of one of its neighbours in the ordering in 
coordinate Cd~i, and it differs from the location of the other neighbour in coordinate Cj where j = ao{r). 
Thus, j = ao(r) is the unique j < d—1 such that Cj{f) ^ Cj{r). The entrance gate of subregion S{r) is at 
T(r)((0, 0, 0)), which is the point p = {po, ...,pd-i) such that pi = ^{ci{r) + mi{r)) for < i < d; the 
exit gate is at r(r)((l, 0, 0)), which is the point q — {qo, qd-i) such that qj = ^{cj{r) + l~mj{r)) for 
j = ao(r), and = ^{ci{r)+mi{r)) — pi for all i G {0, d — 1}\ {ao(r)}. We have to show (i) pi — Xi{r) 
and (ii) = Xi{f) for all z G {0, d— 1}. For i < d— 1, we can rewrite (i) as: Ci{r)+mi(r) = Ci(f)+Ci(r), 
and therefore: mi{r) — Ci{f ). This corresponds directly to the definition of mi{r). For i < d— 1 and i ^ j, 
we can rewrite (ii) as: Ci(r) + mi(r) = Ci(r) + Ci{r) = Ci{r) + Ci(r); again, this corresponds directly to the 
definition of TOi(r). For j = ao(r), we rewrite (ii) as: Cj{r) + l — mj{r) = Cj{r) + Cj{r) = Cj{r) + 1 — Cj{f), 
and therefore: mj{r) — Cj{f). Again, this corresponds directly to the definition of mi{r). Finally, for 
i = d — 1, we rewrite both (i) and (ii) as: Ci{r) + mi(r) = 1. This, too, corresponds directly to the 
definition of mi{r), QED. □ 

Lemma 9 The Butz-Moore curves and the harmonious Hilbert curves defined above are standard Hilbert 
curves. 

Proof: What we have to prove is that the permutations a in the definition of the Butz-Moore curves and 
the harmonious Hilbert curves satisfy the definition of do in standard Hilbert curves. 

For the Butz-Moore curves this is trivial, as their permutations are defined in terms of a given aq. 

For the harmonious Hilbert curves the proof is more involved. First, when r = or r = (l)d, all bits 
of r are equal. Thus, by the definition of harmonious Hilbert curves, ad-i{r) — {d — 1) — {d — 1) = 0, 
which is indeed consistent with the definition of standard Hilbert curves. Second, consider the case when 
< r < 2*^ — 1 and r is odd, and let i have the lowest value such that ri ^ fi. Then i < d — 1, = 0, 
and r^+i, ...,rd^i = 1. By the definition of harmonious Hilbert curves, we will have ai{r) = 0, and by 
the definition of standard Hilbert curves, we will have Ci{r) ^ Ci{r), and thus, indeed, ai(r) — 0. Third, 
consider the case when < r < 2^^ — 1 and r is even, and let i have the lowest value such that ^ fi. 
Then i < d — 1, = 1, and r^+i, ...,rd^i — 0. By the definition of harmonious Hilbert curves, we will 
have ai{r) — 0, and by the definition of standard Hilbert curves, we will have Ci{r) 7^ Ci(f), and thus, 
indeed, ai{r) = 0. □ 

Corollary 1 The Butz-Moore curves and the harmonious Hilbert curves defined above are mono-Hilbert 
curves. 



5.4 Proof of interdimensional consistency 

In this section we will prove that the harmonious Hilbert curves are interdimensionally consistent. Let 



fd denote the d-dimensional harmonious Hilbert curve. As observed in Section 3.2 to prove that the 



harmonious Hilbert curves are interdimensionally consistent, it suffices to show that each curve fd (for 
d > 1) shows fd~i on each front face. In fact, we will prove something stronger, namely that each curve 
fd shows fd-i on each front face, and each back face Fl with < i < d — 1 shows fd~i mirrored in 
dimension i (note that this dimension corresponds to dimension i + 1 in the d-dimensional space) . This 
leaves F^_^ as the only face for which we do not prove anything about the visible order. 
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Below, we write / for fd, and by c, a, to, o, M , and r we denote the location, permutation, reflection, 
translation, transformation matrix and transformation functions that define /. We write /' for fd-i, and 
by c', a', m', o', M' and r' we denote the location, permutation and reflection, translation, transformation 
matrix and transformation functions that define /'. The proof that each curve / shows /' (sometimes 
mirrored, as specified above) on all faces except F^_^ goes by induction on increasing level of refinement. 
More precisely, for a level of refinement £, we consider the recursive balanced subdivision of the (d— 1)- 
dimensional unit hypercube U' into 2^'-''"^' regions, and prove the following induction hypothesis: 

(i) the visible order of / on any face F' — Ff with Q < i < d — 1 visits these regions in the same order 
as /'; 

(ii) the visible order of / on any face F' = with Q < i < d~ 1 visits these regions in the same order 
as /' mirrored in dimension i. 

The proof follows the same general approach as the proof of consistency for the sets of mono-Wunderlich 
curves discussed in Section |4l 



Base case 1=1 As a base case, consider refinement level £ = 1, and a face F' = F^ . Recall that 
/ is defined based on subdivision of a d-dimensional unit hypercube — which we will denote U — into 2*^ 
regions S'(O), S{{l)d)- This induces a subdivision of F' into 2^*^^ regions i^'(O), each of 
which is a (d — 1) -dimensional face of one of the regions S'(O), ^((l)^)- Let F'(0), be 
indexed in the order in which the c?-dimensional regions that contain them appear in S{0), S{{l)d)- 
Thus, there is a monotonously increasing function a : {0, {l)d-i] ^ {0, such that F'{r') is a 
face of S{a{r')). Note that for r' G {0, (1)^-1}, the locations of the regions F'ir') within U' are given 
by X c(iT(r')) = X 2RGC''((T(r')). Furthermore, the co-domain of a consists of exactly those values r such 
that Ci{r) = 2RGcf (r) = k. Thus, the locations of the regions F'{r') within [/', in order of increasing r', 
form the sequence '^'^^RGC'^, which, by Lemma |6j is exactly the sequence 2rgc'^~^ if /c = 0, and the 
sequence q2rgc''~^ if fc = 1 and i < d — 1. Hence, the visible order of / on any face F' = F^ visits 
the regions of refinement level £ = 1 in the same order as /', or /' mirrored in dimension i, respectively, 
thus establishing the base case for our induction. 



Induction For the inductive step, we have to be more careful than with the mono-Wunderlich curves. 
Before, all we had to do is proving that the transformations that map the curve / through U to the 
section of / through S{a{r')), also transform the order /' shown on a face F*{r') into the section of 
/' that should be visible on F'Jr'). More concretely, we showed that we had to prove that we have 
X m{a{r')) = m! ir') (Equation [l]) and * o.ifi'''')) — oli^') (Equation [2|. However, with the harmonious 
Hilbert curves, not all faces show /': some faces show a mirrored version, and one face (i^J_i) docs not 
show /' in any way. Denote the face of U that maps to F'ir') by F*{r'). We will use r as a shorthand 
for o'(r'). Recall that we now have: 

F*(r'] ^ 

Therefore, we have to show: 

• that we never have k ^ mi{r) while ai{r) = d — 1 (never use the 'broken' face); 

• that the reflections TO(r) undo any mirroring of the order on F*{r') and apply any necessary 
mirroring on F'; 

• that the permutations match: x '^(r) = a'{r') (this is no different from the situation with mono- 



Wunderlich curves in Section 4.4) 



Recall that, by the definition of F' as F^ , we have Ci{r) = k. 

We will first establish that we have r' = ^r. If fc = 0, then, by definition of r = o'(r'), we have 
c(r) = *=''c'(r'), and hence, 2RGcf{r)= and * 2RGC'*(r) = 2RGC''~^(r'). By Lemma IS] we now 
have ri = and thus, by Lemma 6 part (i), case (a), we have r' = ^r. Similarly, if fc = 1 and 

i < d — 1, then, by definition of r = a\r'), we have c(r) = Qc'(r'), and hence, 2RGcf(r) = 1 and 
*2RGC'*(r) = (^2RGC'^~i(r'). By Lemma [S] we now have ^ and thus, by Lemma 6 part (i), 

case (c), we have r' = ^r. 



29 



Never the broken face For the sake of contradiction, suppose k ^ mi{r) and ai{r) = d — 1. We 
distinguish three cases: r = 0; r > and i < d — 1; and i = d — 1. 

First, consider the case of r = 0. Then k = Ci{r) and mi{r) = Ci(r) by definition, so we cannot have 
k 7^ mi{r). 

Second, consider the case of r > and i < d — 1. By the definition of the permutations a{r) in 
harmonious Hilbert curves, it follows from ai{r) — d~l that = r^-i. However, by the definition of the 
reflections mi{r) in standard Hilbert curves, we also have Ci{f) — mi{r). Since Ci{r) ^ k =/= jni{r) = Ci{f), 
this implies (by Lemmajs]) that i is the smallest index such that ^ fi and therefore r.^ = 1 and rd-i = 0. 
This contradicts the observation that = r^-i-, so we cannot have k ^ mi(r) and ai{r) — d — I. 

Third, consider the case of i = — 1. In this case, our induction hypothesis only concerns the case 
k = 0, and thus we have Cd-i{r) = 0, ad-i{r) = d ^ 1 and md-i{r) = 1. By the definition of the 
permutations a(r) in harmonious Hilbert curves, ad-i{r) ^ d—l if and only if r.i ^ r^-i for alH < rf— 1. 
But then, by Lemmajs] Cd-i(r) = \rd^i — rd-2 \ = 1, which contradicts our observation that Cd-i{r) = 0. 
So in this case we cannot have k ^ mi{r) and ai{r) — d — \ either. 

Reflections We need to prove that the reflections m{r) undo any mirroring of the order on F* {r') and 
apply any necessary mirroring on F' . Let j be a.i(r) and let g be aj+i(r), that is, g is defined such that 
ag{r) — aj(r) + l. We define r7i(r) as the reflections of S'(r) with application (or equivalently, cancellation) 
of the mirroring of the order of /' as shown on F*{r'). Note that the order of /' as shown on F*{r') is 
mirrored only if mi(r) ^ k, and if so, it is mirrored in coordinate aj(r), which is coordinate ai{r) + 1 
in d-dimensional space on F*{r'), and this coordinate is mapped to coordinate g on F'. Therefore, we 
have m(r) = m{r) if mi(r) = k, and m(r) = Qm(r) if rriiir) ^ k. Furthermore, we define m! {r') 
as the reflections of F'(r) with the necessary mirroring of F' , that is, m'{r') — m!{r') if A: = 0, and 
rh'ir') = Qm'{r') if fc = 1. What we need to prove is now: ^ 'rh{r) = rh'{r'). 

We deflne c as follows: if mi(r) = fc, then c = c; if rai{r) ^ k, then c= qC. Thus, m(r) is expressed 
in terms of c(r) and c(f) just as m(r) is expressed in terms of c(r) and c(f) in the definition of standard 
Hilbert curves. Similarly, we define c' as follows: if fc = 0, then c' = c'; if fc = 1, then c' — qc'. Thus, 
m! (r') is expressed in terms of c'(r') and c'(r') just as m(r) is expressed in terms of c(r) and c(f ) in the 
definition of standard Hilbert curves. 

Note that & gives the locations of the subregions of U' that touch F', in the order in which these 
subregions are visited by /; thus we have, for any q' £ {0, {l)d-i}: 

c'{q')= lc{a{q')). (5) 

We are now ready to prove l^rh{r) = m! {r'). We distinguish four cases: 

(1) mi{r) — k and r = 0; 

(2) mi{r) — k and r > and i < d — 1; 

(3) mi{r) — k and r > and i = d — 1; 

(4) mi{r) 7^ k and r' = 0; 

(5) mi{r) 7^ k and r' > 0. 

In case (1), we have m{r) — m{r) (because mi{r) — fc), we have Ci{r) = (because r — 0), and thus, 
k = and m!{r') — m'{r'). Since r' — — we have m{r) = m{r) = c{r) — and m!{r') = m'{r') — 
c'(r') — 0, and therefore ^ 'Ti(r) = rh'{r'), QED. 

In case (2), we have m(r) = m{r) (because mi(r) — fc). We have Ci{f) — mi{r) = fc. Thus, the 
subregion preceding S{r) in the d-dimensional order also contributes the region preceding F'{r') on F' , 
of which the location is given by c'(r') (note that if fc = 1, then /' as visible on F' is mirrored in 
dimension i). Therefore we have c'(r') = ^ c(f) and c'(r') = * c(r). Therefore, for z < i we have: 
^2:(''') = c'j,(f') = Cz(f) = TO2(r); for i < z < d — 2 we have: m'^{r') — c'^{f') — c^-|_i(r) — mz+i{r); and 
finally we have fT?4-2(^') = 1 ^ "^d-2(''') ~ ^ ~ Cd-ii^) — fnd-i{r). Thus, m'{r') = * m(r) — ^ m{r), 
QED. 

In case (3), by the definition of the permutations m in a standard Hilbert curve and the definition 
of the case mi{r) = fc, we have Cd~i{r) = 1 — 'md-i{r) = 1 — fc, but this contradicts Ci[r) — fc, so this 
cannot occur. 
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In case (4), we have m,;(r) ^ k, and therefore ?fi(r) = Qm{r), where g is such that ag{r) — ai{r) + 1. 
Since r' — 0, we have that r equals the smallest q such that Ci{q) = k. If fc = 0, this implies r = 0, but 
then, by the definition of reflections in a standard Hilbert curve, we would have rriiir) = Ci(r) = = fc, 
contradicting the definition of case (4). Therefore, k = 1, and the induction hypothesis only concerns 
cases with i < d — 1. By Lemmajsj we now have r = (0)il0(0)(i_i_2 and r = (0)i01(l)d-i-2- Applying 
the definitions of the permutations a and the reflections in m in standard and harmonious Hilbert curves, 
we find ai{r) = 0, .g = ai(r) = d - 1, c(r) = (0)ill(0)d-i-2, c(f) = (0)i01(0)d_j_2; if i < d - 2 we have 
m{r) — (0)i01(0)ci_i_3l; if ? = d — 2 we have m(r) = (0)^00; therefore we have m(r) — '^Q^m(r) = 
{0)M{0)d-^-2 in all cases, and ^ m(r) = (0)a(0)d-»-2 - ©(O)d-i = m'(0), QED. 

In case (5), we have mi{r) ^ fc, and therefore m(r) = Qm{r), where g is such that ag{r) — ai{r) + 1. 
Since r' > 0, we also have r > 0, and both r and f' exist. We have Ci{f) = m.i{r) ^ k. Thus, the 
subregion preceding S{r) in the d-dimensional order does not contain the region preceding F'{r') on F' . 
Let p' be f' and let p be (j{p')- Note that we have Ci{p) = k, Ci{p), Ci{f) — 1 ~ k, Ci{r) = k, and: 



. c(f) 



c(r). 



(6) 



By Lemma [3j as the rank increases from p to r, digit i must be the most significant digit that changes 
from to 1 between p and p, and the next time it changes from to 1 is between f and r. Therefore, 
between p and r, there is a unique more significant digit h that changes from to 1, and it does so 
exactly once. Thus we get the following values for the digits of p, p, f and r: 



Xo,---,Xh-l Xh Xh+l, ■■■^Xi-i 



p anything 10 1 

p same as p 1 1 

f same as p 1 1 

r same as p 1 1 

Note that we may have h — i ~ 1 or i ^ d ~ 1, so digits Xh+i, and digits Xi+i, ...,Xd-i might 

not exist. When rank digit h changes as the rank increases from p to r, the corresponding digit in the 
location changes as well. Thus, h must have the unique value such that c'f^ip') c'^(r'), and we get: 



Since Ci{p) — Ci{r), we have that h also has the unique value such that c/i(p) =/= Ch{r), so: 

c(p) = ec{r) 



(7) 



(8) 



Applying the definition of the permutations a in harmonious Hilbert curves to the numbers in the table. 



we find that we have a^(r) = ai{r) + 1. Thus, since mi{r) = 1 — fc, we have c 



_ h 



(just as r' = * r). Using Equations [sj [6] [t] and [s] and h < i, we get: 

c'(f') = c'(p') = gc'(r') = ^ gc(r) = ^ gc(f) 



Note that p' = ^p 



c(f). 



(9) 



Thus, for z < i we have: rh'^{r') = c^(r') = Cz{f) 
Cz+i{r) = m2+i(r); and finally we have f7'4-2(^') " 
Thus, m'(r') = ^ m(r), QED. 



rhz{r); for i < z < d — 2 we have rh'^\ 



(r') = 1 - Crf_i(r) = 1 - Cd-i(r) = rhd-iir). 



Permutations Recall that we have r' = * r. Therefore, if i < d— 1, we have r^_2 = J^d-i- If « = d— 1, 
the induction hypothesis only concerns the case of A; = 0, that is, Cd-i{r) = 0, and therefore, by Lemmajsj 
we have rd-i = rd-2 = i^d-2- Thus, we always have r'^_2 = ^d-i- 

Recall that the respective inverses of a'(r') and a{r) (where r = cr{r')) can be constructed from the 
identity permutation by moving all indices h such that = '''d-2^ or — Vd-i = t'^_2^ respectively, to 
the front, reversing the order of those indices, and reversing the order of the remaining indices. Inserting 
a digit in r' at position i to obtain r = cr(r'), does not affect which of the other digits are moved to the 
front, since this only depends on the value of r^_2- Thus, inserting a digit in r' at position i to obtain r, 
moving the index i in a(r), and subsequently removing index i from d{r) (which is what ^ a(^) does), does 
not affect the order of the remaining indices relative to each other. Therefore we have a'(r') = * 0'{f')^ 
QED. 
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5.5 Orientation of the Butz-Moore and harmonious Hilbert curves 

By Lemma [2] a d-dimensional mono-Hilbert curve has neutral orientation if and only if there is a 
number £, such that each of the d\ possible permutations of d numbers can be constructed as the 
composition of £ permutations from a. Trivially, each 1-dimensional curve has neutral orientation. Below 
we analyse whether the d-dimensional curves for d > 2 have neutral orientation. 

Butz-Moore curves In the Butz-Moore curves, all permutations are rotations. Only if d = 2, this 
suffices to create all possible permutations by composition. Therefore, only the two-dimensional Butz- 
Moore curve, that is, only the original Hilbert curve, has neutral orientation; the higher-dimensional 
Butz-Moore curves have not. 

Harmonious Hilbert curves In the harmonious Hilbert curves, a((0)(j_2l0) is the permutation that 
swaps the first two elements and then reverses the whole sequence, and a((0)d-i(l)i) is the permutation 
that rotates the sequence i steps to the right and then reverses the whole sequence. This allows us to 
construct any swap of two adjacent elements at positions j — 1 and j by composing the four permutations 
a((0)j+i(l)d-j-i), a((0)d_2lO), a((0)d-j-i(l)j+i) and a((0)rf). Since any permutation can be created 
from any other by at most d(d — l)/2 swaps of adjacent elements, we get that any permutation can be 
created from any other by composing at most 2d{d — 1) permutations from a. However, to satisfy the 
condition in Lemma |2j there should be a depth i such that each permutation can be created from the 
identity permutation by composing exactly £ permutations. In the harmonious Hilbert curves, this can 
indeed be done for £ = 2d{d— 1), since any permutation that is produced from k < d{d— l)/2 swaps can 
be produced by a sequence of Ak permutations producing the swaps, followed by 2d{d — 1) — 4/c reversal 
permutations a((0)d). It follows that the harmonious Hilbert curves have neutral orientation. 

Note that the analysis given above is by no means tight: in reality much less than 2d{d — l)/2 
permutations from a may suffice to produce all possible permutations. For example, in four dimensions 
three permutations suffice, rather than 2d{d — 1) = 24. 

6 First-half- and diagonal-extradimensional space-filling curves 

6.1 Definition of first-half- and diagonal-extradimensionahty 

Recall that a d-dimensional space-filling curve / is extradimensional to a d'-dimensional space-filling 
curve /' with co-domain [/', if the following holds for any valid ordering /' of /' and any choice of ji: 
there is a valid ordering / of / such that for any pair of points a,b G U' we have /'(a) < f'{b) if and only 
if /(liftp(a)) < /(lift^(6)) (Definition ll|. Here lift^ is a function that defines the correspondence between 
points in the co-domains of / and / . the sequence of coordinates of lift^ (p) is simply the sequence of 
coordinates of p, with zeros inserted in certain places according to /i. In this section we will consider two 
variants of extradimensionality: first-half-extradimensionality and diagonal-extradimensionality. 

In the case of first-half-extradimensionality, d = 2c?' and ^ is restricted to the identity function id, 
defined by id{i) — i. Thus, the coordinates of lift^(p) are simply the coordinates of p with d' zeros added 
at the end: 

Definition 3 A 2d' -dimensional space-filling curve f is first-half-extradimensional to a d' -dimensional 
space-filling curve f with co-domain U' , if the following holds for any valid ordering f of f : there is 
a valid ordering f of f such that for any pair of points a, 6 € [/' we have f'{a) < f'{b) if and only if 
7(lift,<j(a)) <7(lift,rf(6)). 

Given a d'-dimensional point p — (po, •■•jPd'-i), we define double (p) as the 2d'-dimensional point 
(go, 92ci'-i) such that qi = j^od d')- In other words, double(p) is obtained from p by concatenating 
two copies of the coordinate sequence of p. For a set of points P, let double(P) be the set IJpeP double(p). 

Definition 4 A 2d' -dimensional space-filling curve f is diagonal-extradimensional to a d'-dimensional 
space-filling curve f with co-domain U' , if the following holds for any valid ordering f of f : there is 
a valid ordering f of f such that for any pair of points a,b d U' we have f'{a) < f'{b) if and only if 
7(double(a)) < 7(double(5)). 
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Inset 4 Building R-trees with extradimensional space-filling curves. 

An R-tree is a balanced tree structure in wiiicii tiie leaves store bounding boxes of d-dimensional objects. Each object is 
stored exactly once, and typically each leaf contains the bounding boxes of a fair number of objects: enough to fill a page on 
disk. As the bounding box of an object, the smallest axis-aligned d-dimensional box that contains the object is used. Each 
leaf or other node of the tree also has a bounding box of its own, which is the smallest axis-aligned d-dimensional box that 
encloses all objects stored in the subtree rooted at that node. In addition, each non-leaf node stores the bounding boxes 
of its children. This structure can be used to answer several types of queries efficiently [5], provided the object bounding 
boxes are distributed over the leaves in such a way that leaves have small bounding boxes. 

One heuristic way to achieve this, is by the use of space-filling curves 1111 . To determine the ordering of two-dimensional 
axis-parallel rectangles (object bounding boxes) in the tree, we represent the rectangles as four-dimensional points, which 
are subsequently sorted by the order of their positions along a four-dimensional space-filling curve. Two ways to represent 
rectangles [aJmin i ^^max] x [j/min i J/max] as four-dimensional points were considered. The xy-representation would use the 
point: 

(^'mim ymin; ^S'niax, ymax) , 

whereas the cd-representation (centre/dimensions) would use the point: 

(^min "t- ^Jmax 2/min -t- ^max 
2 ' 2 

Since in practice, many rectangles are small, xy-representation leads to many points where the first and the third coordinate 
are almost identical, and the second and the fourth coordinate are almost identical. The other solution, cd-representation, 
leads to many points whose last two coordinates are almost zero. In previous work we argued and showed evidence that it 
is important that such four-dimensional points are ordered as much as possible as if the represented rectangles were simply 
ordered according to the positions of their centre points along a good two-dimensional space-filling curve |9|. This can 
be achieved by using four-dimensional curves that are diagonal- or first-half-extradimensional to Hilbert's two-dimensional 
curve. 

The results in this paper could be used to try this approach to building R-trees in higher dimensions (for example, using 
six-dimensional curves for three-dimensional boxes), or to experiment with 3-regular curves instead of 2-regular curves for 
this application. 



In previous work, we found that certain applications call for space-filling curves that are first-half- 
and/or diagonal-extradimensional to a given space- filling curve with good locality-preserving proper- 
ties [9]. For details, see Inset [4] 

Any 2d'-dimensional Peano curve, coil curve, half-coil curve, Meurthe curve, or harmonious Hilbert 
curve (see Sections[4]and[5| is first-half-extradimensional, but generally not diagonal-extradimensional, to 
the d'-dimensional Peano curve, coil curve, half-coil curve, Meurthe curve, or harmonious Hilbert curve, 
respectively. In this section we will show how, given any d'-dimensional ^-regular order-preserving vertex- 
continuous mono-curve /' with b £ {2,3} and /'(O) = (0,0,. ..,0), one can construct a 2c?'-dimensional 
6-regular vertex-continuous mono-curve / that is both first-half-extradimensional and diagonal-extra- 
dimensional to /'. This may provide a useful alternative to the curves from the previous sections 
if diagonal-extradimensionality is desired, or if one wants to base the 2(i'-dimensional curve on a d'- 
dimensional curve that is different from those presented in the previous sections. 

6.2 Monotone space-filling curves 

As a building block in our construction of first-half- and diagonal-extradimensional curves, we will use 
two-dimensional curves that have certain special properties. We will denote these curves by /i2 and h^, 
where /12 is the two-dimensional Hilbert curve (Figure [lOj left), and is the 5I-curve (Figure 
Lemma [TO] states that /i2 and have the properties which we will need to prove the correctness of our 
construction in Section [6.31 

Definition 5 An ordering f of a space-filling curve f is monotone from a to b, if f{p) increases mono- 
tonically as a point p moves in a straight line from a to b. 

Lemma 10 The curves h2 and have valid orderings h2 and h^ that are monotone from (0, 0) to (1, 0) 
and from (0,0) to (1,1). 

Proof; We will prove the lemma for the Hilbert curve /12; the proof for the H-curve h^ is very similar 
and is left as an exercise for the interested reader. In fact, we will prove a stronger claim for ft,2, namely 
that h2 has a valid ordering /12 that is monotone from (0, 0) to (0, 1), from (0, 1) to (1, 1), from (1, 1) to 
(1,0), from (0,0) to (1,0), and from (0,0) to (1,1). 
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Figure 10: Hilbert's curve (left) and the 5I-eurve (right), defined graphically and by means of a table. 



We will prove the lemma by induction on increasing level of refinement of the subdivision of the unit 
square into subsquares. More precisely, for a level of refinement £, we consider the recursive balanced 
subdivision of the unit square into 4^ square regions, and prove the following induction hypothesis: there 
is a valid ordering h2 that: 

(1) puts the regions touching the left edge of the unit square in order from bottom (0,0) to top (0, 1); 

(2) puts the regions touching the top edge of the unit square in order from left (0, 1) to right (1, 1); 

(3) puts the regions touching the right edge of the unit square in order from top (1, 1) to bottom (1, 0); 

(4) puts the regions touching the bottom edge of the unit square in order from left (0, 0) to right (1,0); 

(5) puts the regions on the ascending diagonal of the unit square in order from bottom (0, 0) to top 
(1,1)- 

As a base case, we take £ = 1. For this base case, the induction hypothesis can easily be verified by 
inspection of Figure [lO) 

We will now prove that if the induction hypothesis holds for refinement level £~ I, then it also holds 
for refinement level £. Note that, as far as the left edge, the bottom edge and the diagonal are concerned, 
the induction hypothesis is invariant under the transformation applied in the first (lower left) quadrant 
of and, as far as the right and the bottom edge are concerned, the induction hypothesis is invariant 
under the transformation applied in the last (lower right) quadrant of /i2. Thus h2 places: 

• the regions touching the lower half of the left edge of the unit square before the regions touching 
the upper half of the left edge of the unit square (by the induction hypothesis for level 1); 

• the regions touching the lower half of the left edge of the unit square in order from bottom to top 
(by the induction hypothesis for level £ — 1); 

• the regions touching the upper half of the left edge of the unit square in order from bottom to top 
(by the induction hypothesis for level £ — 1). 

As a result, /i2 puts all level-^ regions touching the left edge of the unit square in order from bottom to 
top. In a similar way, one can verify the order induced by ^2 on the regions touching the top, right, and 
bottom edges and on the ascending diagonal of the unit square. □ 



6.3 Constructing first-half- and diagonal-extradimensional curves 

To be able to present and analyse our construction, we need some notation. When / : [0, 1] — > [0, 1]'* is a 
d-dimensional space-filling curve, let /^(t) be coordinate i of f{t), where coordinates are numbered from 
to d-1. In other words, f{t) = {f°\t),f^\t), /'''"^'(i))- For two d'-dimensional points p and q, we 
let p\q be the 2ci'-dimensional point obtained by concatenating the coordinate sequences of p and q, that 
is, p\q = {po, ...,pd'-i,qo, ■■■,qd'-i)- For a 2c?'-dimensional point p, we let denote the c?'-dimensional 
point obtained by taking the first d' coordinates of p, that is, — {po, ...,pd'-i)- Similary, p^ is the 
point obtained by taking the last d' coordinates of p, that is, p^ = {pd', ■■■,P2d'-i)- We use not only 
to denote the number zero, but also to denote the origin of the coordinate system, in any number of 
dimensions. 
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Now, given a d'-dimensional &-regular space-filling curve /' with b £ {2,3}, we can define a 2d'- 
dimcnsional space-filling curve / by: 

/w = /'(4°'w)i/'(4"w) 

To put it differently: for k e {0, 1} and i e {0, ...,d' ~ 1}, we define as (^f' (*))■ With a 

slight abuse of terminology and notation, we can think of / as the composition of /' and /i;,, denoted 
/' o h\,. Below we will prove the following: 

Theorem 3 ///' is an order-preserving vertex- continuous mono-curve and /'(O) = 0, then f = f ° hb 
is a b-regular vertex-continuous space-filling curve, and f is first-half- extradimensional and diagonal- 
extradimensional to /'. //, in addition, f is symmetric, then f is an order-preserving mono-curve. 

In Section [6.4| we will see how to implement a comparison operator for /. We will now prove Theorem[3] 



in four parts (Lemmas 11 ■ 14) 



Lemma 11 // /' is a order-preserving b-regular vertex- continuous mono-curve, then f = f ° is a 
b-regular vertex- continuous space-filling curve. 

Proof: For each point p € [0,1]^^^ there are values a,b £ [0,1] such that /'(a) = and f'{b) = p'~ 
(because /' is a surjective function from [0, 1] to [0, 1]'' ), and there is a value t such that h2{t) = (a, b) 
(because h2 is a surjective function from [0, 1] to [0, 1]^). Therefore / is a surjective function from [0, 1] 
to [0, l]^'' . Furthermore, since /' and hb are continuous, so is /. 

It remains to prove that / is 6-regular. Consider level d'£ of refinement of hj,. At this level of 
refinement, hi, maps each interval [r/b"^'^ ^, (r + l)/6^'^ for r G {0, b^"^ ^ — 1}, to a square of width 
1/b'^^. The horizontal range of the square, that is, the possible values of that are attained in this 
square, correspond exactly to the pre- image of one of the b'^ ^ level-^ regions of /', thus defining the region 
occupied in 2(i'-dimensional space with respect to the first d' dimensions. Similarly, the vertical range of 
the square corresponds exactly to the pre-image of a region of /', thus defining the region occupied in 
2d'-dimensional space with respect to the last d' dimensions. Thus, each interval [r/h'^'^ ^, (r-f- \)/b'^'^ ^] is 
mapped exactly to the Cartesian product of two d'-dimensional hypercube regions of width 1/6^, which 
is a 2d'-dimensional hypercube region of width 1/b^. This implies the hierarchical structure of a 5-regular 
curve. □ 

Lemma 12 // /' is a symmetric b-regular mono-curve, then f — f o hj, is an order-preserving b-regular 
mono- curve. 

Proof: Extending our notation from Section [2. 3( we denote the subregions, locations, permutations, and 
reflections of a curve / by 5-^, , and , respectively. Let hb-e be the level i of refinement of hb, 
at which level regions are numbered from to 6^^ — 1, and the region coordinates c^*"^ are between 
and 6^ — 1 (inclusive). The transformation of / in region r g {0, ...,6^^^ — 1} is the result of the 
successive effects of (i) the permutation a'*f' <''(r), (ii) the reflections m'*'' '*' (r), (iii) the permutations 
(c^' (r)) for k e {0, 1}, and (iv) the reflections (Cfc' {r)). Note that there are no reversals, since 
/' is symmetric, and therefore, order-preserving. 

(i) The permutation a'^*- (r) is either 01 or 10: the latter results in a rotation of / that consists in 
swapping coordinates 0, d' — 1 with coordinates d' , 2d' — 1. 

(ii) Because /' is symmetric, there is a transformation p (consisting of rotation and/or reflections) 
such that /'(I — h^b^{t)) = p{f'{h^^\t))). A reflection mQ*"*' (r) results in replacing the coordinates 

= fihfHt)) by /'(I - hf\t)) = pifihf\t)) = pirit)) for ^ G {0, d - 1}, which results 
in a transformation of / that consists in applying p to the first d' coordinates while leaving the last 
d' coordinates unaffected; thus, this transformation of / consists of rotations and reflections in 2d'- 
dimensional space. Similarly, a reflection rn^* (r) results in a transformation of / that consists in 
applying p to the last d' coordinates while leaving the first d' coordinates unaffected. 

(iiijiv) Similar to p, the permutations a-^ (c)!'' '*' (t")) and the reflections (c^*" (j')) result in trans- 
formations of / that consist of rotations and reflections. 

Thus, the transformation of / in region r e {0, b^'^ — 1} consists of rotations and reflections, and 
/ is an order-preserving mono-curve. □ 
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Table 3: The four-dimensional curve /12 o created from the composition of two-dimensional Hilbert 
curves. This curve is first-half- and diagonal-extradimensional to the two-dimensional Hilbert curve. 



rank 


loc. 


perm. 


rcfl. 


exit 


rank 


loc. 


perm. 


rcfl. 


exit 


0000 


0000 


1032 


0000 


4(0,1,0,0) 
1(0,1,0,1) 
1(0,1,1,1) 

5(0,0,1,1) 


1000 


1111 


2301 


0000 


4(1,1,2,1) 
1(2,1,2,1) 
5(2,0,2,1) 
1(2,0,1,1) 


0001 


0100 


2310 


0000 


1001 


1110 


0132 


0011 


0010 


0101 


2301 


0000 


1010 


1010 


1032 


1111 


0011 


0001 


1023 


0110 


1011 


1011 


3201 


1010 


0100 


0011 


3201 


0000 


5(0,0,2,1) 

1(0,1,2,1) 
|(1,1,2,1) 
5(1,1,1,1) 


1100 


1001 


1023 


1010 


f(2, 1,1,1) 
1(2,1,0, 1) 
5(2,1,0,0) 
5(2,0,0,0) 


0101 


0010 


1032 


0011 


1101 


1101 


2301 


1010 


0110 


0110 


0132 


0011 


1110 


1100 


2310 


1001 


0111 


0111 


2301 


1010 


1111 


1000 


1032 


1100 



Note that if /' is not symmetric, then / is not necessarily a mono-curve. A reflection of /if, in 
coordinate k e {0, 1} results in replacing (/if ' (t)) by fl''d'+i](i _ h^^\t)) for i e {0, ...,d' - 1}. In 

general, this transformation cannot be expressed as a simple permutation, reflection, or reversal of the 
d-dimensional curve /. For example, if 6 = 3 and /' = /13 (the 5I-curve of Figure 10), then / is not a 
mono-curve. The same might be the case if /' is a Meurthe curve (see Section |4]) . 

However, if /' is symmetric (for example, a Butz-Moore Hilbert curve or a harmonious Hilbert curve 
from Section [5] or a Peano curve, a coil curve, or a half-coil curve from Section |4]), then / will be an 
order-preserving mono- cu rve. For example, if /' is a standard Hilbert curve, then (using the notation 

from the proof of Lemma 12) a reflection mj!'' (r) leads to a reflection in axis Sq {cl"-"' (r)). Thus we 
get: 



a 



f 

kd' 

f 

kd' + 
/ 

''kd' + 



m 



(r) 
(r) 



- J' 



= af(4-^-'(r)) + d'. 
= m{ (cj!"-'''(r)) 
= l-mf(4-'"(r)) 



if ml^"' (r) = or af (cj!^""^' (r)) 7^ 
if ml^"' (r) 1 and af (Cfc'"' (^)) = 



Table [s] shows the resulting curve definition if /' is the two-dimensional Hilbert curve, and thus, the 
resulting curve / is /12 o /12. Note that / is not a standard Hilbert curve as defined in Section 5.1 since 
the locations of the subregions do not conform to the binary reflected Gray code. 

Lemma 13 The curve f is first-half-extradimensional to /'. 

Proof; Let /' be any valid ordering of /. By Definition [sj we have to prove that there is a valid ordering 
/ of / such that for any pair of points a,b £ U' we have /'(a) < /'(6) if and only if /(a|0) < f{b\0). 

Let hb be a valid ordering of /ib which is monotone from (0,0) to (1,0) (such an ordering exists by 
Lemma 10). We now define / by f{p) — hi,{f'{p^), f'{p^)). Recall that /' is a 6-regular curve with 
the entrance gate in the origin, and therefore /'(O) = 0. Thus, we have /(a|0) = /ib((/'(a), /'(O))) — 
h{(fia),0)) and 7(6J0) = M_(f (6), f (0))) = 7^((f (&), 0))^. By the monotonicity of 7^ from (0,0) to 
(1, 0), it follows that /(a|0) < /(&|0) if and only if /'(a) < f'{b), QED. □ 

Lemma 14 The curve f is diagonal-extradimensional to /'. 

Proof: Let /' be any valid ordering of /. By Definition [4j we have to prove that there is a valid ordering 
/ of / such that for any pair of points a,b G U' we have /'(a) < /'(6) if and only if f{a\a) < f{b\b). 

Let hhhe a vahd ordering of hb which is monotone from (0,0) to (1, 1) (such an ordering exists by 
Lemma[l0l). We now define / by /(p) = h~b(f'ip'^),TiP^))- Thus, we have J{a\a) = /^((f (a), /'(«))_) and 
f{b\b) = /ifc((/^6), f'{b))). By the monotonicity of h from (0, 0) to (1, 1), it follows that /(a|a) < f{b\b) 
if and only if /'(a) < /'(5), QED. □ 



6.4 Implementing a first-half- and diagonal-extradimensional curve 

In Section [63] we saw how to one can construct a first-half- and diagonal-extradimensional curve /, given 
a /)-regular order-preserving vertex-continuous mono-curve /' that satisfies certain conditions. We also 
saw that if /' is symmetric, then / is an order-preserving mono-curve, which could be implemented using 
Algorithm [2?2| 
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We will now see another method to implement a comparison operator for /, which can also be 
applied if /' is not symmetric and which may also be more efficient. We illustrate this method by 
implementing first-half- and diagonal-extradimensional curves based on the Butz-Moore generalization 
of Hilbert curves. Our goal is to compare points according to a valid ordering / of /. Considering the 
construction of /, we can implement this by computing f{p) as /i2(/'(p^), f'{p^))- To do so, we rewrite 
Algorithm |5 . 1 1 as a comparator object for Butz-Moore curves that is initialized with the coordinates of 
two d'-dimensional points p' and g', and provides a method to compute f'{p') and f'{q') one digit a time. 
A comparison operator for / is then implemented by initializing two comparator objects, one with 
and , and one with p'~ and g'~, and then simply computing /12, extracting digits from the comparator 
objects rather than directly from the coordinates of p and q. 

The details are as follows. Each object has member variables forward, rotation, reflected, i, iMinusOne, 
newRotation, pLocal and qLocal. To initialize an object, we run Algorithm |6.1[ calling it with two d'- 



dimensional points p and q; to extract a rank digit (one for p and one for q), we run Algorithm 6.2 The 
extraction method works correctly up to and including the first time the rank digit differs between p 
and q (digits extracted afterwards may be wrong). The algorithm for the 2(i'-dimensional curve is given 
as Algorithm |6.3| 



Algorithm 6.1: Initialization of a comparator object for a Butz-Moore curve 

1 pLocal -k— p; qLocal g; forward -s— true; rotation -h- 0; for i -s— to d' — 1 do reflected[i] false 

2 i ^ rotation; iMinusOne -s— d' — 1: newRotation iMinusOne 



Algorithm 6.2: Extracting a single rank digit from a comparator object for a Butz-Moore curve 

1 pDigit ExTRACTFlRSTDlGlT(pLoca/[i]); qDigit ^ ExTRACTFlRSTDlGlT(gLoca/[i]) 

2 if reflected[i] then 

3 1^ pDigit -S— (1 — pDigit); qDigit ^ (1 — qDigit) 

4 if forward then 

5 1^ pRank 4— pDigit; qRank <~ qDigit 

6 else 

7 1^ pRank -S— (1 — pDigit); qRank -S— (1 — qDigit) 

8 if pDigit = 1 then 

9 1^ forward 4— -^forward; reflected[i] 4— -^reflected[i]; newRotation <— iMinusOne 

10 iMinusOne 4— i; « 4— (i + 1) mod d' 

11 if i = rotation then 

12 reflected[iMinusOne] 4— -^reflected[iMinusOne] 

13 reflected[newRotation] 4— ^reflected[newRotation] 

14 rotation 4— newRotation 

15 z 4— rotation; iMinusOne 4— (z — 1) mod d'; newRotation 4— iMinusOne 

16 return [pRank, qRank) 



7 Remaining questions for theory and practice 

Building R-trees As explained in Inset |4j the curves in this paper could be used to build R-trees 
in higher dimensions (for example, using six-dimensional curves for three-dimensional boxes), or to 
experiment with 3-regular curves instead of 2-regular curves for this application. Some first, rough 
experiments, indicate that it may depend a lot on the distribution of the data whether the new curves 
offer any advantages over previously known curves — but it seems unlikely that the new curves have 
serious disadvantages in this application. Further experiments would need to be done. 

Fast implementations To get most out of the space-filling curves presented in this paper, one may 
want to design highly-optimized implementations of comparison operators. As compared to the 3- 
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Algorithm 6.3: Algorithm to compare the position of two points p and q along a 2(i'-dimensional 
curve that is first-half- and diagonal-extradimcnsional to a d'-dimensional Butz-Moore curve. 

1 initialize a comparator object compare[0] with points and 

2 initialize a comparator object compare[l] with points p^ and 

3 forward ^ true; rotation ^ 0; for i to 1 do reflected[i] ^ false 

4 repeat 

5 
6 
7 
8 
9 



10 
11 
12 

13 
14 
15 
16 
17 



(i — 1) mod 2; newRotation ^ iMinusOne 



i <— rotation; iMinusOne 
repeat 

extract (pDigit, qDigit) from compare[i] 
if reflected[i] then 
1^ pDigit -s— (1 — pDigit); qDigit -S— (1 — qDigit) 

if pDigit < qDigit then return forward; else if pDigit > qDigit then return - 
if pDigit — 1 then 

1^ forward -s— -^forward; reflected[i\ -s— ^re/Zecte(i[i]; newRotation -s— iMinusOne 

iMinusOne -s— i; « -s— (i + 1) mod 2 
until i = rotation 

reflected[iMinusOne] ^re/!ecterf[iMmMsOrae] 
reflected[newRotation\ <— -^reflected[newRotation] 
rotation ^ newRotation 

18 until a/Z remaining digit strings p[0], ...,p[2d' — 1] and ^[O], (7[2(i' — 1] are empty 

19 return false 



^forward 



regular curves, the 2-regular curves have the advantage that many calculations can be implemented 
by means of bit shifting and/or exclusive-or operations on many bits at once. A possible advantage of 
interdimensionally consistent curves, is that a highly-optimized (possibly hard- wired) implementation 
of a d-dimensional curve, for some fixed value of d, can also be used as an implementation of each d'- 
dimensional curve with d' < d. If the cost of working with more dimensions than necessary is mitigated 
by parallel circuitry, then, being able to use a single hard-wired d-dimensional implementation for any 
d' < d, may reduce overhead. 

Proving that the harmonious Hilbert curves are interdimensionally consistent At this point, 
I would like to pose the main result of this paper as an open question. In Section |5.4[ I proved that 
the harmonious Hilbert curves consitute an interdimensionally consistent set. However, I believe the 
proof is hard to understand and to verify. In particular, I do not find the current proof more convincing 
than several earlier attempts at a proof, which I only dismissed as erroneous because they led to the 
conclusion that the harmonious Hilbert curves could not be interdimensionally consistent. Indeed, these 
earlier attempts had to be wrong, or my exhaustive search algorithm to find the curves and the algorithms 
I used to test the implementations all had to be wrong. I implemented comparison operators for all curves 
from Section |4] and [5] for d e {2,3,4,5} and verified interdimensional consistency with respect to the 
ordering of subregions down to at least the third level of refinement. Although now, everything seems to 
fall into place and I am convinced that my implementations and Section [5. 4| do not contain any mistakes 
that cannot be repaired, it would still be good if we could find a simpler proof. 

Clearing up interdimensional consistency under tie-breaking Tie-breaking is a subtle issue 
which was touched on in Sections [T] and [2] but not resolved completely. Space-filling curves are surjective, 
but not bijective — in fact, they cannot be bijective. Often this issue is dealt with by considering the 
level-^ approximation of the space-filling curve — say, the Hilbert curve /12 , for concreteness — which is the 
polygonal curve whose vertices are the lower left corners of the 4^ successive squares at the £-th level of 
refinement. Clearly, for any £, the level-£ approximation visits each point of the unit square at most once, 
and in the limit for £ — >■ 00, the approximation visits all points of the unit square. Unfortunately, this 
does not imply that each point is visited only once: this approach offers no way around the fact that, for 
example, /i2(l/6) = h2{l/2). This can be seen by writing out /i2(l/6) and /i2(l/2) in binary notation: 
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/i2(l/110) = /i2(0.00101010...) = (0.0111..., 0.0111...) = (0.1,0.1), and /i2(l/10) = /i2(0.1) = (0.1,0.1). 
This makes the inverse of h2 ih-defined. 

In the definitions and proofs in this paper, we deah with the ambiguity as follows. We allow any 
'inverse' /i2 such that for all points p in the unit hypercube, we have /i2(^2(p)) = P', and we prove that no 
matter which inverse /12 we choose, we can choose at least one valid inverse for each higher-dimensional 
harmonious Hilbert curve so that we get an interdimensionally consistent set. 

However, our implementations of comparison operators are essentially based on level-£ approxima- 
tions. The idea is that this implements the tie-breaking method that always assigns each point p to 
the quadrant that lies to its upper right. In fact, this approach would not work correctly if we would 
accept coordinates that end with an infinite sequence 111...: it would put the points (0.0111,0.0111), 
(0, 1) and (0.1,0.1) in this order, even though the first and the last point are the same. Our algorithms 
do work correctly as long as the coordinates of each point p are represented as binary numbers with a 
finite number of digits. Recall that our algorithms do not accept the number 1 for a coordinate. 

The above observations imply that the proofs of interdimensional consistency for the curves as such, 
are not necessarily valid for the curves with the specific choice of inverses as implemented in our algo- 
rithms. The fact that our implementations of mono-Wunderlich curves are interdimensionally consistent 



could still be established easily by inspecting the algorithms (see Section 4.2). For harmonious Hilbert 
curves the situation is less clear. The algorithm for harmonious Hilbert curves is much harder to analyse. 
The proof that harmonious Hilbert curves are interdimensionally consistent relies on the proper handling 
of points with coordinates equal to 1, but the implementation does not accept the number 1 and has not 
been verified to work correctly when given the alternative representation 0.111... of the same number. 
Although everything seems to indicate that this is ultimately just a matter of presentation, it is one that 
would be good to resolve. 



Interdimensionally consistent sets of mono-Hilbert curves with interior origins? In Sec- 
tion |5.2| I had to give up on extending the harmonious Hilbert curves to filling the full d-dimensional 
space (rather than only the space of points with non-negative coordinates) while retaining interdimen- 
sional consistency. The reason was that I could not find a way to move the origin into the interior of 
the d-dimensional unit hypercube, while maintaining that the c?-dimensional harmonious Hilbert curve 
shows the (d — l)-dimensional curve on the intersection of the unit hypercube with any axis-parallel 
{d — l)-dimensional plane through the origin. I conjecture that this is, in fact, impossible. One could 
also try to find out if an interdimensionally consistent set of 2-regular vertex-continuous space-filling 
curves — not necessarily mono-curves — could exist. 



How unique are sets of interdimensionally consistent curves? The fact that the three-, four-, 
five-, six- and seven-dimensional harmonious Hilbert curves emerged as the unique results of an auto- 
mated search, leads me to conjecture that the harmonious Hilbert curves actually constitute the only 
interdimensionally consistent set of mono-Hilbert curves. The catch is that, for d-dimensional curves 
with d > 4, the search was restricted to symmetric curves in which the subregions are traversed in the 
order of binary refiected Gray codes. Although this seems to be the most natural choice, and although 
this is indeed the only solution for d — 2 and d 3, I cannot be sure that no alternative solutions for 
d > 4 were ruled out by these restrictions. 

While I believe that the mono-Hilbert curves may contain only one interdimensionally consistent set, 
we already saw that the mono-Wunderlich curves contain at least four interdimensionally consistent sets. 
There must be many more, if only because the subregion in the centre can be rotated freely without 
affecting the curves shown on the faces of the unit hypercube. One could try to find a full characterization 
of all interdimensionally consistent sets of mono-Wunderlich curves — possibly with the added condition 
that the back faces also show the (d — l)-dimensional curves. 
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